*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2012, by the GROMACS development team, led by
+ * Copyright (c) 2012,2014, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
* 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 "nb_generic.h"
-#include "nrnb.h"
-#include "nonbonded.h"
-#include "nb_kernel.h"
+#include <math.h>
+#include "gromacs/gmxlib/nonbonded/nb_kernel.h"
+#include "gromacs/legacyheaders/nonbonded.h"
+#include "gromacs/legacyheaders/nrnb.h"
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/legacyheaders/types/simple.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/utility/fatalerror.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;
int ewitab;
real ewtabscale, eweps, sh_ewald, ewrt, ewtabhalfspace;
real * ewtab;
- real rcoulomb2, rvdw, rvdw2, sh_invrc6;
+ real rcoulomb2, rvdw, rvdw2, sh_dispersion, sh_repulsion;
real rcutoff, rcutoff2;
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];
rcoulomb2 = fr->rcoulomb*fr->rcoulomb;
rvdw = fr->rvdw;
rvdw2 = rvdw*rvdw;
- sh_invrc6 = fr->ic->sh_invrc6;
+ 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 */
fvdw = (vvdw_rep-vvdw_disp)*rinvsq;
if (fr->vdw_modifier == eintmodPOTSHIFT)
{
- vvdw = (vvdw_rep-c12*sh_invrc6*sh_invrc6)*(1.0/12.0)-(vvdw_disp-c6*sh_invrc6)*(1.0/6.0);
+ vvdw = (vvdw_rep + c12*sh_repulsion)/12.0 - (vvdw_disp + c6*sh_dispersion)/6.0;
}
else
{
fvdw = (br*vvdw_rep-vvdw_disp)*rinvsq;
if (fr->vdw_modifier == eintmodPOTSHIFT)
{
- vvdw = (vvdw_rep-cexp1*exp(-cexp2*rvdw))-(vvdw_disp-c6*sh_invrc6)/6.0;
+ vvdw = (vvdw_rep-cexp1*exp(-cexp2*rvdw))-(vvdw_disp + c6*sh_dispersion)/6.0;
}
else
{
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 */
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff