Implemented nbnxn LJ switch functions

[alexxy/gromacs.git] / src / gromacs / mdlib / nbnxn_kernels / nbnxn_kernel_ref_inner.h

diff --git a/src/gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_ref_inner.h b/src/gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_ref_inner.h
index a4adc825b29d6c5b590a01fa54d9c33d0afd25a4..6c9c89dec45c10843330b984ed1142bf319bf03a 100644 (file)
--- a/src/gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_ref_inner.h
+++ b/src/gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_ref_inner.h
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2012,2013, by the GROMACS development team, led by
+ * Copyright (c) 2012,2013,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.
@@ -69,7 +69,11 @@
             real rsq, rinv;
             real rinvsq, rinvsix;
             real c6, c12;
-            real FrLJ6 = 0, FrLJ12 = 0, VLJ = 0;
+            real FrLJ6 = 0, FrLJ12 = 0, frLJ = 0, VLJ = 0;
+#if defined VDW_FORCE_SWITCH || defined VDW_POT_SWITCH
+            real r, rsw;
+#endif
+
 #ifdef CALC_COULOMB
             real qq;
             real fcoul;
@@ -90,9 +94,7 @@
              * (if appropriate) the i and j indices are
              * unsuitable for this kind of inner loop. */
             real skipmask;
-#ifdef VDW_CUTOFF_CHECK
-            real skipmask_rvdw;
-#endif
+
 #ifdef CHECK_EXCLS
             /* A multiply mask used to zero an interaction
              * when that interaction should be excluded
@@ -151,26 +153,71 @@
             {
                 rinvsix = interact*rinvsq*rinvsq*rinvsq;
 
-#ifdef VDW_CUTOFF_CHECK
-                skipmask_rvdw = (rsq < rvdw2);
-                rinvsix      *= skipmask_rvdw;
-#endif
-
                 c6      = nbfp[type_i_off+type[aj]*2  ];
                 c12     = nbfp[type_i_off+type[aj]*2+1];
                 FrLJ6   = c6*rinvsix;
                 FrLJ12  = c12*rinvsix*rinvsix;
-                /* 6 flops for r^-2 + LJ force */
-#ifdef CALC_ENERGIES
-                VLJ     = (FrLJ12 - c12*sh_invrc6*sh_invrc6)/12 -
-                    (FrLJ6 - c6*sh_invrc6)/6;
+                frLJ    = FrLJ12 - FrLJ6;
+                /* 7 flops for r^-2 + LJ force */
+#if defined CALC_ENERGIES || defined VDW_POT_SWITCH
+                VLJ     = (FrLJ12 + c12*ic->repulsion_shift.cpot)/12 -
+                    (FrLJ6 + c6*ic->dispersion_shift.cpot)/6;
+                /* 7 flops for LJ energy */
+#endif
+
+#if defined VDW_FORCE_SWITCH || defined VDW_POT_SWITCH
+                /* Force or potential switching from ic->rvdw_switch */
+                r       = rsq*rinv;
+                rsw     = r - ic->rvdw_switch;
+                rsw     = (rsw >= 0.0 ? rsw : 0.0);
+#endif
+#ifdef VDW_FORCE_SWITCH
+                frLJ   +=
+                    -c6*(ic->dispersion_shift.c2 + ic->dispersion_shift.c3*rsw)*rsw*rsw*r
+                    + c12*(ic->repulsion_shift.c2 + ic->repulsion_shift.c3*rsw)*rsw*rsw*r;
+#if defined CALC_ENERGIES
+                VLJ    +=
+                    -c6*(-ic->dispersion_shift.c2/3 - ic->dispersion_shift.c3/4*rsw)*rsw*rsw*rsw
+                    + c12*(-ic->repulsion_shift.c2/3 - ic->repulsion_shift.c3/4*rsw)*rsw*rsw*rsw;
+#endif
+#endif
+
+#if defined CALC_ENERGIES || defined VDW_POT_SWITCH
+                /* Masking shoule be done after force switching,
+                 * but before potential switching.
+                 */
                 /* Need to zero the interaction if r >= rcut
                  * or there should be exclusion. */
                 VLJ     = VLJ * skipmask * interact;
-                /* 9 flops for LJ energy */
+                /* 2 more flops for LJ energy */
+#endif
+
+#ifdef VDW_POT_SWITCH
+                {
+                    real sw, dsw;
+
+                    sw    = 1.0 + (swV3 + (swV4+ swV5*rsw)*rsw)*rsw*rsw*rsw;
+                    dsw   = (swF2 + (swF3 + swF4*rsw)*rsw)*rsw*rsw;
+
+                    frLJ  = frLJ*sw - r*VLJ*dsw;
+                    VLJ  *= sw;
+                }
+#endif
+
 #ifdef VDW_CUTOFF_CHECK
-                VLJ    *= skipmask_rvdw;
+                /* Mask for VdW cut-off shorter than Coulomb cut-off */
+                {
+                    real skipmask_rvdw;
+
+                    skipmask_rvdw = (rsq < rvdw2);
+                    frLJ         *= skipmask_rvdw;
+#ifdef CALC_ENERGIES
+                    VLJ    *= skipmask_rvdw;
 #endif
+                }
+#endif
+
+#ifdef CALC_ENERGIES
 #ifdef ENERGY_GROUPS
                 Vvdw[egp_sh_i[i]+((egp_cj>>(nbat->neg_2log*j)) & egp_mask)] += VLJ;
 #else
@@ -244,8 +291,8 @@
             if (i < UNROLLI/2)
 #endif
             {
-                fscal = (FrLJ12 - FrLJ6)*rinvsq + fcoul;
-                /* 3 flops for scalar LJ+Coulomb force */
+                fscal = frLJ*rinvsq + fcoul;
+                /* 2 flops for scalar LJ+Coulomb force */
             }
 #ifdef HALF_LJ
             else
@@ -254,7 +301,7 @@
             }
 #endif
 #else
-            fscal = (FrLJ12 - FrLJ6)*rinvsq;
+            fscal = frLJ*rinvsq;
 #endif
             fx = fscal*dx;
             fy = fscal*dy;