Sort all includes in src/gromacs

[alexxy/gromacs.git] / src / gromacs / gmxlib / nonbonded / nb_generic.c

diff --git a/src/gromacs/gmxlib/nonbonded/nb_generic.c b/src/gromacs/gmxlib/nonbonded/nb_generic.c
index e81ca49723b92d10986691e2e6c6cf9f50d8bdef..55fdb28a01f334790ca3605d4b49a544fbe2d7cf 100644 (file)
--- a/src/gromacs/gmxlib/nonbonded/nb_generic.c
+++ b/src/gromacs/gmxlib/nonbonded/nb_generic.c
@@ -34,21 +34,19 @@
  * 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,
@@ -76,10 +74,10 @@ 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;
@@ -97,6 +95,7 @@ gmx_nb_generic_kernel(t_nblist *                nlist,
     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];
@@ -120,6 +119,11 @@ gmx_nb_generic_kernel(t_nblist *                nlist,
     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)
     {
@@ -154,7 +158,7 @@ gmx_nb_generic_kernel(t_nblist *                nlist,
 
     bExactElecCutoff    = (fr->coulomb_modifier != eintmodNONE) || fr->eeltype == eelRF_ZERO;
     bExactVdwCutoff     = (fr->vdw_modifier != eintmodNONE);
-    bExactCutoff        = bExactElecCutoff || bExactVdwCutoff;
+    bExactCutoff        = bExactElecCutoff && bExactVdwCutoff;
 
     if (bExactCutoff)
     {
@@ -172,7 +176,7 @@ gmx_nb_generic_kernel(t_nblist *                nlist,
     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;
 
@@ -182,6 +186,7 @@ gmx_nb_generic_kernel(t_nblist *                nlist,
     shiftvec            = fr->shift_vec[0];
     vdwparam            = fr->nbfp;
     ntype               = fr->ntype;
+    vdwgridparam        = fr->ljpme_c6grid;
 
     for (n = 0; (n < nlist->nri); n++)
     {
@@ -222,7 +227,7 @@ gmx_nb_generic_kernel(t_nblist *                nlist,
             velec            = 0;
             vvdw             = 0;
 
-            if (bExactCutoff && rsq > rcutoff2)
+            if (bExactCutoff && rsq >= rcutoff2)
             {
                 continue;
             }
@@ -251,6 +256,10 @@ gmx_nb_generic_kernel(t_nblist *                nlist,
                         /* 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:
@@ -309,8 +318,8 @@ gmx_nb_generic_kernel(t_nblist *                nlist,
                 }
                 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 */
@@ -391,6 +400,29 @@ gmx_nb_generic_kernel(t_nblist *                nlist,
                         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;
@@ -408,8 +440,8 @@ gmx_nb_generic_kernel(t_nblist *                nlist,
                 }
                 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 */