* To help us fund GROMACS development, we humbly ask that you cite
* the research papers on the package. Check out http://www.gromacs.org.
*/
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
#include <math.h>
-#include "types/simple.h"
-#include "vec.h"
-#include "typedefs.h"
+#include "gromacs/legacyheaders/types/simple.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/legacyheaders/typedefs.h"
#include "nb_generic.h"
-#include "nrnb.h"
+#include "gromacs/legacyheaders/nrnb.h"
-#include "nonbonded.h"
-#include "nb_kernel.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/legacyheaders/nonbonded.h"
+#include "nb_kernel.h"
void
gmx_nb_generic_kernel(t_nblist * nlist,
real ix, iy, iz, fix, fiy, fiz;
real jx, jy, jz;
real dx, dy, dz, rsq, rinv;
- real c6, c12, cexp1, cexp2, br;
+ real c6, c12, c6grid, cexp1, cexp2, br;
real * charge;
real * shiftvec;
- real * vdwparam;
+ real * vdwparam, *vdwgridparam;
int * shift;
int * type;
real * fshift;
real rswitch_elec, rswitch_vdw, 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;
+ real ewclj, ewclj2, ewclj6, ewcljrsq, poly, exponent, sh_lj_ewald;
gmx_bool bExactElecCutoff, bExactVdwCutoff, bExactCutoff;
x = xx[0];
rvdw2 = rvdw*rvdw;
sh_dispersion = fr->ic->dispersion_shift.cpot;
sh_repulsion = fr->ic->repulsion_shift.cpot;
+ sh_lj_ewald = fr->ic->sh_lj_ewald;
+
+ ewclj = fr->ewaldcoeff_lj;
+ ewclj2 = ewclj*ewclj;
+ ewclj6 = ewclj2*ewclj2*ewclj2;
if (fr->coulomb_modifier == eintmodPOTSWITCH)
{
bExactElecCutoff = (fr->coulomb_modifier != eintmodNONE) || fr->eeltype == eelRF_ZERO;
bExactVdwCutoff = (fr->vdw_modifier != eintmodNONE);
- bExactCutoff = bExactElecCutoff || bExactVdwCutoff;
+ bExactCutoff = bExactElecCutoff && bExactVdwCutoff;
if (bExactCutoff)
{
eps = 0.0;
eps2 = 0.0;
- /* 3 VdW parameters for buckingham, otherwise 2 */
+ /* 3 VdW parameters for Buckingham, otherwise 2 */
nvdwparam = (ivdw == GMX_NBKERNEL_VDW_BUCKINGHAM) ? 3 : 2;
table_nelements = 12;
shiftvec = fr->shift_vec[0];
vdwparam = fr->nbfp;
ntype = fr->ntype;
+ vdwgridparam = fr->ljpme_c6grid;
for (n = 0; (n < nlist->nri); n++)
{
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 */
vvdw = vvdw_disp + vvdw_rep;
break;
+
+ case GMX_NBKERNEL_VDW_LJEWALD:
+ /* LJ-PME */
+ rinvsix = rinvsq*rinvsq*rinvsq;
+ ewcljrsq = ewclj2*rsq;
+ exponent = exp(-ewcljrsq);
+ poly = exponent*(1.0 + ewcljrsq + ewcljrsq*ewcljrsq*0.5);
+ c6 = vdwparam[tj];
+ c12 = vdwparam[tj+1];
+ c6grid = vdwgridparam[tj];
+ vvdw_disp = (c6-c6grid*(1.0-poly))*rinvsix;
+ vvdw_rep = c12*rinvsix*rinvsix;
+ fvdw = (vvdw_rep - vvdw_disp - c6grid*(1.0/6.0)*exponent*ewclj6)*rinvsq;
+ if (fr->vdw_modifier == eintmodPOTSHIFT)
+ {
+ vvdw = (vvdw_rep + c12*sh_repulsion)/12.0 - (vvdw_disp + c6*sh_dispersion - c6grid*sh_lj_ewald)/6.0;
+ }
+ else
+ {
+ vvdw = vvdw_rep/12.0-vvdw_disp/6.0;
+ }
+ break;
+
default:
gmx_fatal(FARGS, "Death & horror! No generic VdW interaction for ivdw=%d.\n", ivdw);
break;
}
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 */