Sort all includes in src/gromacs
[alexxy/gromacs.git] / src / gromacs / mdlib / nbnxn_kernels / simd_2xnn / nbnxn_kernel_simd_2xnn_outer.h
index eb7409906c7f13bf02cccf98ebc54826fa49757d..21a63ab4aa7b6fcdd899101711d8c75690025dd6 100644 (file)
@@ -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.
@@ -75,7 +75,7 @@
     real               *nbfp_ptr;
     int                 n, ci, ci_sh;
     int                 ish, ish3;
-    gmx_bool            do_LJ, half_LJ, do_coul;
+    gmx_bool            do_LJ, half_LJ, do_coul, do_self;
     int                 sci, scix, sciy, sciz, sci2;
     int                 cjind0, cjind1, cjind;
     int                 ip, jp;
     real       *vctp[UNROLLI];
 #endif
 
-    gmx_mm_pr  shX_S;
-    gmx_mm_pr  shY_S;
-    gmx_mm_pr  shZ_S;
-    gmx_mm_pr  ix_S0, iy_S0, iz_S0;
-    gmx_mm_pr  ix_S2, iy_S2, iz_S2;
-    gmx_mm_pr  fix_S0, fiy_S0, fiz_S0;
-    gmx_mm_pr  fix_S2, fiy_S2, fiz_S2;
+    gmx_simd_real_t  shX_S;
+    gmx_simd_real_t  shY_S;
+    gmx_simd_real_t  shZ_S;
+    gmx_simd_real_t  ix_S0, iy_S0, iz_S0;
+    gmx_simd_real_t  ix_S2, iy_S2, iz_S2;
+    gmx_simd_real_t  fix_S0, fiy_S0, fiz_S0;
+    gmx_simd_real_t  fix_S2, fiy_S2, fiz_S2;
     /* We use an i-force SIMD register width of 4 */
-    /* The pr4 stuff is defined in nbnxn_kernel_simd_utils.h */
-    gmx_mm_pr4 fix_S, fiy_S, fiz_S;
+    /* The simd4 stuff might be defined in nbnxn_kernel_simd_utils.h */
+    gmx_simd4_real_t fix_S, fiy_S, fiz_S;
 
-    gmx_mm_pr  diagonal_jmi_S;
+    gmx_simd_real_t  diagonal_jmi_S;
 #if UNROLLI == UNROLLJ
-    gmx_mm_pb  diagonal_mask_S0, diagonal_mask_S2;
+    gmx_simd_bool_t  diagonal_mask_S0, diagonal_mask_S2;
 #else
-    gmx_mm_pb  diagonal_mask0_S0, diagonal_mask0_S2;
-    gmx_mm_pb  diagonal_mask1_S0, diagonal_mask1_S2;
+    gmx_simd_bool_t  diagonal_mask0_S0, diagonal_mask0_S2;
+    gmx_simd_bool_t  diagonal_mask1_S0, diagonal_mask1_S2;
 #endif
 
-    unsigned      *exclusion_filter;
-    gmx_exclfilter filter_S0, filter_S2;
+    unsigned            *exclusion_filter;
+    gmx_exclfilter       filter_S0, filter_S2;
 
-    gmx_mm_pr      zero_S = gmx_set1_pr(0);
+    gmx_simd_real_t      zero_S = gmx_simd_set1_r(0.0);
 
-    gmx_mm_pr      one_S = gmx_set1_pr(1.0);
-    gmx_mm_pr      iq_S0 = gmx_setzero_pr();
-    gmx_mm_pr      iq_S2 = gmx_setzero_pr();
-    gmx_mm_pr      mrc_3_S;
+    gmx_simd_real_t      one_S = gmx_simd_set1_r(1.0);
+    gmx_simd_real_t      iq_S0 = gmx_simd_setzero_r();
+    gmx_simd_real_t      iq_S2 = gmx_simd_setzero_r();
+
+#ifdef CALC_COUL_RF
+    gmx_simd_real_t      mrc_3_S;
 #ifdef CALC_ENERGIES
-    gmx_mm_pr      hrc_3_S, moh_rc_S;
+    gmx_simd_real_t      hrc_3_S, moh_rc_S;
+#endif
 #endif
 
 #ifdef CALC_COUL_TAB
     /* Coulomb table variables */
-    gmx_mm_pr   invtsp_S;
-    const real *tab_coul_F;
+    gmx_simd_real_t   invtsp_S;
+    const real       *tab_coul_F;
 #ifndef TAB_FDV0
-    const real *tab_coul_V;
+    const real       *tab_coul_V;
 #endif
-    int         ti0_array[2*GMX_SIMD_WIDTH_HERE], *ti0;
-    int         ti2_array[2*GMX_SIMD_WIDTH_HERE], *ti2;
+    /* Thread-local working buffers for force and potential lookups */
+    int               ti0_array[2*GMX_SIMD_REAL_WIDTH], *ti0 = NULL;
+    int               ti2_array[2*GMX_SIMD_REAL_WIDTH], *ti2 = NULL;
 #ifdef CALC_ENERGIES
-    gmx_mm_pr   mhalfsp_S;
+    gmx_simd_real_t   mhalfsp_S;
 #endif
 #endif
 
 #ifdef CALC_COUL_EWALD
-    gmx_mm_pr beta2_S, beta_S;
+    gmx_simd_real_t beta2_S, beta_S;
 #endif
 
 #if defined CALC_ENERGIES && (defined CALC_COUL_EWALD || defined CALC_COUL_TAB)
-    gmx_mm_pr  sh_ewald_S;
+    gmx_simd_real_t  sh_ewald_S;
+#endif
+
+#if defined LJ_CUT && defined CALC_ENERGIES
+    gmx_simd_real_t   p6_cpot_S, p12_cpot_S;
+#endif
+#ifdef LJ_POT_SWITCH
+    gmx_simd_real_t   rswitch_S;
+    gmx_simd_real_t   swV3_S, swV4_S, swV5_S;
+    gmx_simd_real_t   swF2_S, swF3_S, swF4_S;
+#endif
+#ifdef LJ_FORCE_SWITCH
+    gmx_simd_real_t   rswitch_S;
+    gmx_simd_real_t   p6_fc2_S, p6_fc3_S;
+    gmx_simd_real_t   p12_fc2_S, p12_fc3_S;
+#ifdef CALC_ENERGIES
+    gmx_simd_real_t   p6_vc3_S, p6_vc4_S;
+    gmx_simd_real_t   p12_vc3_S, p12_vc4_S;
+    gmx_simd_real_t   p6_6cpot_S, p12_12cpot_S;
+#endif
+#endif
+#ifdef LJ_EWALD_GEOM
+    real              lj_ewaldcoeff2, lj_ewaldcoeff6_6;
+    gmx_simd_real_t   mone_S, half_S, lje_c2_S, lje_c6_6_S, lje_vc_S;
 #endif
 
 #ifdef LJ_COMB_LB
-    const real *ljc;
+    const real       *ljc;
 
-    gmx_mm_pr   hsig_i_S0, seps_i_S0;
-    gmx_mm_pr   hsig_i_S2, seps_i_S2;
+    gmx_simd_real_t   hsig_i_S0, seps_i_S0;
+    gmx_simd_real_t   hsig_i_S2, seps_i_S2;
 #else
 #ifdef FIX_LJ_C
-    real        pvdw_array[2*UNROLLI*UNROLLJ+GMX_SIMD_WIDTH_HERE];
-    real       *pvdw_c6, *pvdw_c12;
-    gmx_mm_pr   c6_S0, c12_S0;
-    gmx_mm_pr   c6_S2, c12_S2;
+    real              pvdw_array[2*UNROLLI*UNROLLJ+GMX_SIMD_REAL_WIDTH];
+    real             *pvdw_c6, *pvdw_c12;
+    gmx_simd_real_t   c6_S0, c12_S0;
+    gmx_simd_real_t   c6_S2, c12_S2;
 #endif
 
-#ifdef LJ_COMB_GEOM
-    const real *ljc;
+#if defined LJ_COMB_GEOM || defined LJ_EWALD_GEOM
+    const real       *ljc;
 
-    gmx_mm_pr   c6s_S0, c12s_S0;
-    gmx_mm_pr   c6s_S1, c12s_S1;
-    gmx_mm_pr   c6s_S2 = gmx_setzero_pr(), c12s_S2 = gmx_setzero_pr();
-    gmx_mm_pr   c6s_S3 = gmx_setzero_pr(), c12s_S3 = gmx_setzero_pr();
+    gmx_simd_real_t   c6s_S0, c12s_S0;
+    gmx_simd_real_t   c6s_S2  = gmx_simd_setzero_r();
+    gmx_simd_real_t   c12s_S2 = gmx_simd_setzero_r();
 #endif
 #endif /* LJ_COMB_LB */
 
-    gmx_mm_pr  vctot_S, Vvdwtot_S;
-    gmx_mm_pr  sixth_S, twelveth_S;
+    gmx_simd_real_t  vctot_S, Vvdwtot_S;
+    gmx_simd_real_t  sixth_S, twelveth_S;
 
-    gmx_mm_pr  avoid_sing_S;
-    gmx_mm_pr  rc2_S;
+    gmx_simd_real_t  avoid_sing_S;
+    gmx_simd_real_t  rc2_S;
 #ifdef VDW_CUTOFF_CHECK
-    gmx_mm_pr  rcvdw2_S;
-#endif
-
-#ifdef CALC_ENERGIES
-    gmx_mm_pr  sh_invrc6_S, sh_invrc12_S;
-
-    /* cppcheck-suppress unassignedVariable */
-    real       tmpsum_array[2*GMX_SIMD_WIDTH_HERE], *tmpsum;
-#endif
-#ifdef CALC_SHIFTFORCES
-    /* cppcheck-suppress unassignedVariable */
-    real       shf_array[2*GMX_SIMD_WIDTH_HERE], *shf;
+    gmx_simd_real_t  rcvdw2_S;
 #endif
 
     int ninner;
     int npair = 0;
 #endif
 
-#if defined LJ_COMB_GEOM || defined LJ_COMB_LB
+#if defined LJ_COMB_GEOM || defined LJ_COMB_LB || defined LJ_EWALD_GEOM
     ljc = nbat->lj_comb;
-#else
+#endif
+#if !(defined LJ_COMB_GEOM || defined LJ_COMB_LB)
     /* No combination rule used */
     nbfp_ptr    = (4 == nbfp_stride) ? nbat->nbfp_s4 : nbat->nbfp;
 #endif
 
     /* Load j-i for the first i */
-    diagonal_jmi_S    = gmx_load_pr(nbat->simd_2xnn_diagonal_j_minus_i);
+    diagonal_jmi_S    = gmx_simd_load_r(nbat->simd_2xnn_diagonal_j_minus_i);
     /* Generate all the diagonal masks as comparison results */
 #if UNROLLI == UNROLLJ
-    diagonal_mask_S0  = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
-    diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
-    diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
-    diagonal_mask_S2  = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
+    diagonal_mask_S0  = gmx_simd_cmplt_r(zero_S, diagonal_jmi_S);
+    diagonal_jmi_S    = gmx_simd_sub_r(diagonal_jmi_S, one_S);
+    diagonal_jmi_S    = gmx_simd_sub_r(diagonal_jmi_S, one_S);
+    diagonal_mask_S2  = gmx_simd_cmplt_r(zero_S, diagonal_jmi_S);
 #else
 #if 2*UNROLLI == UNROLLJ
-    diagonal_mask0_S0 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
-    diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
-    diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
-    diagonal_mask0_S2 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
-    diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
-    diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
-    diagonal_mask1_S0 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
-    diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
-    diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
-    diagonal_mask1_S2 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
+    diagonal_mask0_S0 = gmx_simd_cmplt_r(zero_S, diagonal_jmi_S);
+    diagonal_jmi_S    = gmx_simd_sub_r(diagonal_jmi_S, one_S);
+    diagonal_jmi_S    = gmx_simd_sub_r(diagonal_jmi_S, one_S);
+    diagonal_mask0_S2 = gmx_simd_cmplt_r(zero_S, diagonal_jmi_S);
+    diagonal_jmi_S    = gmx_simd_sub_r(diagonal_jmi_S, one_S);
+    diagonal_jmi_S    = gmx_simd_sub_r(diagonal_jmi_S, one_S);
+    diagonal_mask1_S0 = gmx_simd_cmplt_r(zero_S, diagonal_jmi_S);
+    diagonal_jmi_S    = gmx_simd_sub_r(diagonal_jmi_S, one_S);
+    diagonal_jmi_S    = gmx_simd_sub_r(diagonal_jmi_S, one_S);
+    diagonal_mask1_S2 = gmx_simd_cmplt_r(zero_S, diagonal_jmi_S);
 #endif
 #endif
 
         exclusion_filter = nbat->simd_exclusion_filter2;
     }
 
-    /* Here we cast the exclusion masks from unsigned * to int * or
-     * real *.  Since we only check bits, the actual value they
-     * represent does not matter, as long as both mask and exclusion
-     * info are treated the same way.
+    /* Here we cast the exclusion filters from unsigned * to int * or real *.
+     * Since we only check bits, the actual value they represent does not
+     * matter, as long as both filter and mask data are treated the same way.
      */
     filter_S0 = gmx_load_exclusion_filter(exclusion_filter + 0*2*UNROLLJ*filter_stride);
     filter_S2 = gmx_load_exclusion_filter(exclusion_filter + 1*2*UNROLLJ*filter_stride);
 
+#ifdef CALC_COUL_RF
+    /* Reaction-field constants */
+    mrc_3_S  = gmx_simd_set1_r(-2*ic->k_rf);
+#ifdef CALC_ENERGIES
+    hrc_3_S  = gmx_simd_set1_r(ic->k_rf);
+    moh_rc_S = gmx_simd_set1_r(-ic->c_rf);
+#endif
+#endif
+
 #ifdef CALC_COUL_TAB
     /* Generate aligned table index pointers */
     ti0 = prepare_table_load_buffer(ti0_array);
     ti2 = prepare_table_load_buffer(ti2_array);
 
-    invtsp_S  = gmx_set1_pr(ic->tabq_scale);
+    invtsp_S  = gmx_simd_set1_r(ic->tabq_scale);
 #ifdef CALC_ENERGIES
-    mhalfsp_S = gmx_set1_pr(-0.5/ic->tabq_scale);
+    mhalfsp_S = gmx_simd_set1_r(-0.5/ic->tabq_scale);
 #endif
 
 #ifdef TAB_FDV0
 #endif /* CALC_COUL_TAB */
 
 #ifdef CALC_COUL_EWALD
-    beta2_S = gmx_set1_pr(ic->ewaldcoeff_q*ic->ewaldcoeff_q);
-    beta_S  = gmx_set1_pr(ic->ewaldcoeff_q);
+    beta2_S = gmx_simd_set1_r(ic->ewaldcoeff_q*ic->ewaldcoeff_q);
+    beta_S  = gmx_simd_set1_r(ic->ewaldcoeff_q);
 #endif
 
 #if (defined CALC_COUL_TAB || defined CALC_COUL_EWALD) && defined CALC_ENERGIES
-    sh_ewald_S = gmx_set1_pr(ic->sh_ewald);
+    sh_ewald_S = gmx_simd_set1_r(ic->sh_ewald);
+#endif
+
+    /* LJ function constants */
+#if defined CALC_ENERGIES || defined LJ_POT_SWITCH
+    sixth_S      = gmx_simd_set1_r(1.0/6.0);
+    twelveth_S   = gmx_simd_set1_r(1.0/12.0);
+#endif
+
+#if defined LJ_CUT && defined CALC_ENERGIES
+    /* We shift the potential by cpot, which can be zero */
+    p6_cpot_S    = gmx_simd_set1_r(ic->dispersion_shift.cpot);
+    p12_cpot_S   = gmx_simd_set1_r(ic->repulsion_shift.cpot);
+#endif
+#ifdef LJ_POT_SWITCH
+    rswitch_S = gmx_simd_set1_r(ic->rvdw_switch);
+    swV3_S    = gmx_simd_set1_r(ic->vdw_switch.c3);
+    swV4_S    = gmx_simd_set1_r(ic->vdw_switch.c4);
+    swV5_S    = gmx_simd_set1_r(ic->vdw_switch.c5);
+    swF2_S    = gmx_simd_set1_r(3*ic->vdw_switch.c3);
+    swF3_S    = gmx_simd_set1_r(4*ic->vdw_switch.c4);
+    swF4_S    = gmx_simd_set1_r(5*ic->vdw_switch.c5);
+#endif
+#ifdef LJ_FORCE_SWITCH
+    rswitch_S = gmx_simd_set1_r(ic->rvdw_switch);
+    p6_fc2_S  = gmx_simd_set1_r(ic->dispersion_shift.c2);
+    p6_fc3_S  = gmx_simd_set1_r(ic->dispersion_shift.c3);
+    p12_fc2_S = gmx_simd_set1_r(ic->repulsion_shift.c2);
+    p12_fc3_S = gmx_simd_set1_r(ic->repulsion_shift.c3);
+#ifdef CALC_ENERGIES
+    {
+        gmx_simd_real_t mthird_S  = gmx_simd_set1_r(-1.0/3.0);
+        gmx_simd_real_t mfourth_S = gmx_simd_set1_r(-1.0/4.0);
+
+        p6_vc3_S     = gmx_simd_mul_r(mthird_S,  p6_fc2_S);
+        p6_vc4_S     = gmx_simd_mul_r(mfourth_S, p6_fc3_S);
+        p6_6cpot_S   = gmx_simd_set1_r(ic->dispersion_shift.cpot/6);
+        p12_vc3_S    = gmx_simd_mul_r(mthird_S,  p12_fc2_S);
+        p12_vc4_S    = gmx_simd_mul_r(mfourth_S, p12_fc3_S);
+        p12_12cpot_S = gmx_simd_set1_r(ic->repulsion_shift.cpot/12);
+    }
+#endif
+#endif
+#ifdef LJ_EWALD_GEOM
+    mone_S           = gmx_simd_set1_r(-1.0);
+    half_S           = gmx_simd_set1_r(0.5);
+    lj_ewaldcoeff2   = ic->ewaldcoeff_lj*ic->ewaldcoeff_lj;
+    lj_ewaldcoeff6_6 = lj_ewaldcoeff2*lj_ewaldcoeff2*lj_ewaldcoeff2/6;
+    lje_c2_S         = gmx_simd_set1_r(lj_ewaldcoeff2);
+    lje_c6_6_S       = gmx_simd_set1_r(lj_ewaldcoeff6_6);
+    /* Determine the grid potential at the cut-off */
+    lje_vc_S         = gmx_simd_set1_r(ic->sh_lj_ewald);
 #endif
-
-    q                   = nbat->q;
-    type                = nbat->type;
-    facel               = ic->epsfac;
-    shiftvec            = shift_vec[0];
-    x                   = nbat->x;
-
-    avoid_sing_S = gmx_set1_pr(NBNXN_AVOID_SING_R2_INC);
 
     /* The kernel either supports rcoulomb = rvdw or rcoulomb >= rvdw */
-    rc2_S    = gmx_set1_pr(ic->rcoulomb*ic->rcoulomb);
+    rc2_S    = gmx_simd_set1_r(ic->rcoulomb*ic->rcoulomb);
 #ifdef VDW_CUTOFF_CHECK
-    rcvdw2_S = gmx_set1_pr(ic->rvdw*ic->rvdw);
+    rcvdw2_S = gmx_simd_set1_r(ic->rvdw*ic->rvdw);
 #endif
 
-#ifdef CALC_ENERGIES
-    sixth_S      = gmx_set1_pr(1.0/6.0);
-    twelveth_S   = gmx_set1_pr(1.0/12.0);
-
-    sh_invrc6_S  = gmx_set1_pr(ic->sh_invrc6);
-    sh_invrc12_S = gmx_set1_pr(ic->sh_invrc6*ic->sh_invrc6);
-#endif
+    avoid_sing_S = gmx_simd_set1_r(NBNXN_AVOID_SING_R2_INC);
 
-    mrc_3_S  = gmx_set1_pr(-2*ic->k_rf);
-
-#ifdef CALC_ENERGIES
-    hrc_3_S  = gmx_set1_pr(ic->k_rf);
-
-    moh_rc_S = gmx_set1_pr(-ic->c_rf);
-#endif
-
-#ifdef CALC_ENERGIES
-    tmpsum   = gmx_simd_align_real(tmpsum_array);
-#endif
-#ifdef CALC_SHIFTFORCES
-    shf      = gmx_simd_align_real(shf_array);
-#endif
+    q                   = nbat->q;
+    type                = nbat->type;
+    facel               = ic->epsfac;
+    shiftvec            = shift_vec[0];
+    x                   = nbat->x;
 
 #ifdef FIX_LJ_C
-    pvdw_c6  = gmx_simd_align_real(pvdw_array);
+    pvdw_c6  = gmx_simd_align_r(pvdw_array);
     pvdw_c12 = pvdw_c6 + UNROLLI*UNROLLJ;
 
     for (jp = 0; jp < UNROLLJ; jp++)
         pvdw_c12[2*UNROLLJ+jp] = nbat->nbfp[0*2+1];
         pvdw_c12[3*UNROLLJ+jp] = nbat->nbfp[0*2+1];
     }
-    c6_S0            = gmx_load_pr(pvdw_c6 +0*UNROLLJ);
-    c6_S1            = gmx_load_pr(pvdw_c6 +1*UNROLLJ);
-    c6_S2            = gmx_load_pr(pvdw_c6 +2*UNROLLJ);
-    c6_S3            = gmx_load_pr(pvdw_c6 +3*UNROLLJ);
-
-    c12_S0           = gmx_load_pr(pvdw_c12+0*UNROLLJ);
-    c12_S1           = gmx_load_pr(pvdw_c12+1*UNROLLJ);
-    c12_S2           = gmx_load_pr(pvdw_c12+2*UNROLLJ);
-    c12_S3           = gmx_load_pr(pvdw_c12+3*UNROLLJ);
+    c6_S0            = gmx_simd_load_r(pvdw_c6 +0*UNROLLJ);
+    c6_S1            = gmx_simd_load_r(pvdw_c6 +1*UNROLLJ);
+    c6_S2            = gmx_simd_load_r(pvdw_c6 +2*UNROLLJ);
+    c6_S3            = gmx_simd_load_r(pvdw_c6 +3*UNROLLJ);
+
+    c12_S0           = gmx_simd_load_r(pvdw_c12+0*UNROLLJ);
+    c12_S1           = gmx_simd_load_r(pvdw_c12+1*UNROLLJ);
+    c12_S2           = gmx_simd_load_r(pvdw_c12+2*UNROLLJ);
+    c12_S3           = gmx_simd_load_r(pvdw_c12+3*UNROLLJ);
 #endif /* FIX_LJ_C */
 
 #ifdef ENERGY_GROUPS
         ci               = nbln->ci;
         ci_sh            = (ish == CENTRAL ? ci : -1);
 
-        shX_S = gmx_load1_pr(shiftvec+ish3);
-        shY_S = gmx_load1_pr(shiftvec+ish3+1);
-        shZ_S = gmx_load1_pr(shiftvec+ish3+2);
+        shX_S = gmx_simd_load1_r(shiftvec+ish3);
+        shY_S = gmx_simd_load1_r(shiftvec+ish3+1);
+        shZ_S = gmx_simd_load1_r(shiftvec+ish3+2);
 
 #if UNROLLJ <= 4
         sci              = ci*STRIDE;
         do_LJ   = (nbln->shift & NBNXN_CI_DO_LJ(0));
         do_coul = (nbln->shift & NBNXN_CI_DO_COUL(0));
         half_LJ = ((nbln->shift & NBNXN_CI_HALF_LJ(0)) || !do_LJ) && do_coul;
+#ifdef LJ_EWALD_GEOM
+        do_self = TRUE;
+#else
+        do_self = do_coul;
+#endif
 
 #ifdef ENERGY_GROUPS
         egps_i = nbat->energrp[ci];
             }
         }
 #endif
-#if defined CALC_ENERGIES
+
+#ifdef CALC_ENERGIES
 #if UNROLLJ == 4
-        if (do_coul && l_cj[nbln->cj_ind_start].cj == ci_sh)
+        if (do_self && l_cj[nbln->cj_ind_start].cj == ci_sh)
 #endif
 #if UNROLLJ == 8
-        if (do_coul && l_cj[nbln->cj_ind_start].cj == (ci_sh>>1))
+        if (do_self && l_cj[nbln->cj_ind_start].cj == (ci_sh>>1))
 #endif
         {
-            int  ia;
-            real Vc_sub_self;
+            if (do_coul)
+            {
+                real Vc_sub_self;
+                int  ia;
 
 #ifdef CALC_COUL_RF
-            Vc_sub_self = 0.5*ic->c_rf;
+                Vc_sub_self = 0.5*ic->c_rf;
 #endif
 #ifdef CALC_COUL_TAB
 #ifdef TAB_FDV0
-            Vc_sub_self = 0.5*tab_coul_F[2];
+                Vc_sub_self = 0.5*tab_coul_F[2];
 #else
-            Vc_sub_self = 0.5*tab_coul_V[0];
+                Vc_sub_self = 0.5*tab_coul_V[0];
 #endif
 #endif
 #ifdef CALC_COUL_EWALD
-            /* beta/sqrt(pi) */
-            Vc_sub_self = 0.5*ic->ewaldcoeff_q*M_2_SQRTPI;
+                /* beta/sqrt(pi) */
+                Vc_sub_self = 0.5*ic->ewaldcoeff_q*M_2_SQRTPI;
 #endif
 
-            for (ia = 0; ia < UNROLLI; ia++)
+                for (ia = 0; ia < UNROLLI; ia++)
+                {
+                    real qi;
+
+                    qi = q[sci+ia];
+#ifdef ENERGY_GROUPS
+                    vctp[ia][((egps_i>>(ia*egps_ishift)) & egps_imask)*egps_jstride]
+#else
+                    Vc[0]
+#endif
+                        -= facel*qi*qi*Vc_sub_self;
+                }
+            }
+
+#ifdef LJ_EWALD_GEOM
             {
-                real qi;
+                int  ia;
+
+                for (ia = 0; ia < UNROLLI; ia++)
+                {
+                    real c6_i;
 
-                qi = q[sci+ia];
+                    c6_i = nbat->nbfp[nbat->type[sci+ia]*(nbat->ntype + 1)*2]/6;
 #ifdef ENERGY_GROUPS
-                vctp[ia][((egps_i>>(ia*egps_ishift)) & egps_imask)*egps_jstride]
+                    vvdwtp[ia][((egps_i>>(ia*egps_ishift)) & egps_imask)*egps_jstride]
 #else
-                Vc[0]
+                    Vvdw[0]
 #endif
-                    -= facel*qi*qi*Vc_sub_self;
+                        += 0.5*c6_i*lj_ewaldcoeff6_6;
+                }
             }
+#endif      /* LJ_EWALD */
         }
 #endif
 
         gmx_load1p1_pr(&iy_S2, x+sciy+2);
         gmx_load1p1_pr(&iz_S0, x+sciz);
         gmx_load1p1_pr(&iz_S2, x+sciz+2);
-        ix_S0          = gmx_add_pr(ix_S0, shX_S);
-        ix_S2          = gmx_add_pr(ix_S2, shX_S);
-        iy_S0          = gmx_add_pr(iy_S0, shY_S);
-        iy_S2          = gmx_add_pr(iy_S2, shY_S);
-        iz_S0          = gmx_add_pr(iz_S0, shZ_S);
-        iz_S2          = gmx_add_pr(iz_S2, shZ_S);
+        ix_S0          = gmx_simd_add_r(ix_S0, shX_S);
+        ix_S2          = gmx_simd_add_r(ix_S2, shX_S);
+        iy_S0          = gmx_simd_add_r(iy_S0, shY_S);
+        iy_S2          = gmx_simd_add_r(iy_S2, shY_S);
+        iz_S0          = gmx_simd_add_r(iz_S0, shZ_S);
+        iz_S2          = gmx_simd_add_r(iz_S2, shZ_S);
 
         if (do_coul)
         {
-            gmx_mm_pr facel_S;
+            gmx_simd_real_t facel_S;
 
-            facel_S    = gmx_set1_pr(facel);
+            facel_S    = gmx_simd_set1_r(facel);
 
             gmx_load1p1_pr(&iq_S0, q+sci);
             gmx_load1p1_pr(&iq_S2, q+sci+2);
-            iq_S0      = gmx_mul_pr(facel_S, iq_S0);
-            iq_S2      = gmx_mul_pr(facel_S, iq_S2);
+            iq_S0      = gmx_simd_mul_r(facel_S, iq_S0);
+            iq_S2      = gmx_simd_mul_r(facel_S, iq_S2);
         }
 
 #ifdef LJ_COMB_LB
             nbfp3 = nbfp_ptr + type[sci+3]*nbat->ntype*nbfp_stride;
         }
 #endif
+#endif
+#ifdef LJ_EWALD_GEOM
+        /* We need the geometrically combined C6 for the PME grid correction */
+        gmx_load1p1_pr(&c6s_S0, ljc+sci2+0);
+        if (!half_LJ)
+        {
+            gmx_load1p1_pr(&c6s_S2, ljc+sci2+2);
+        }
 #endif
 
         /* Zero the potential energy for this list */
-        Vvdwtot_S        = gmx_setzero_pr();
-        vctot_S          = gmx_setzero_pr();
+        Vvdwtot_S        = gmx_simd_setzero_r();
+        vctot_S          = gmx_simd_setzero_r();
 
         /* Clear i atom forces */
-        fix_S0           = gmx_setzero_pr();
-        fix_S2           = gmx_setzero_pr();
-        fiy_S0           = gmx_setzero_pr();
-        fiy_S2           = gmx_setzero_pr();
-        fiz_S0           = gmx_setzero_pr();
-        fiz_S2           = gmx_setzero_pr();
+        fix_S0           = gmx_simd_setzero_r();
+        fix_S2           = gmx_simd_setzero_r();
+        fiy_S0           = gmx_simd_setzero_r();
+        fiy_S2           = gmx_simd_setzero_r();
+        fiz_S0           = gmx_simd_setzero_r();
+        fiz_S2           = gmx_simd_setzero_r();
 
         cjind = cjind0;
 
 #define CALC_LJ
         if (half_LJ)
         {
+            /* Coulomb: all i-atoms, LJ: first half i-atoms */
 #define CALC_COULOMB
 #define HALF_LJ
 #define CHECK_EXCLS
             while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
             {
-#include "nbnxn_kernel_simd_2xnn_inner.h"
+#include "gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_inner.h"
                 cjind++;
             }
 #undef CHECK_EXCLS
             for (; (cjind < cjind1); cjind++)
             {
-#include "nbnxn_kernel_simd_2xnn_inner.h"
+#include "gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_inner.h"
             }
 #undef HALF_LJ
 #undef CALC_COULOMB
         }
         else if (do_coul)
         {
+            /* Coulomb: all i-atoms, LJ: all i-atoms */
 #define CALC_COULOMB
 #define CHECK_EXCLS
             while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
             {
-#include "nbnxn_kernel_simd_2xnn_inner.h"
+#include "gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_inner.h"
                 cjind++;
             }
 #undef CHECK_EXCLS
             for (; (cjind < cjind1); cjind++)
             {
-#include "nbnxn_kernel_simd_2xnn_inner.h"
+#include "gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_inner.h"
             }
 #undef CALC_COULOMB
         }
         else
         {
+            /* Coulomb: none, LJ: all i-atoms */
 #define CHECK_EXCLS
             while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
             {
-#include "nbnxn_kernel_simd_2xnn_inner.h"
+#include "gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_inner.h"
                 cjind++;
             }
 #undef CHECK_EXCLS
             for (; (cjind < cjind1); cjind++)
             {
-#include "nbnxn_kernel_simd_2xnn_inner.h"
+#include "gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_inner.h"
             }
         }
 #undef CALC_LJ
 
         /* Add accumulated i-forces to the force array */
         fix_S = gmx_mm_transpose_sum4h_pr(fix_S0, fix_S2);
-        gmx_store_pr4(f+scix, gmx_add_pr4(fix_S, gmx_load_pr4(f+scix)));
+        gmx_simd4_store_r(f+scix, gmx_simd4_add_r(fix_S, gmx_simd4_load_r(f+scix)));
 
         fiy_S = gmx_mm_transpose_sum4h_pr(fiy_S0, fiy_S2);
-        gmx_store_pr4(f+sciy, gmx_add_pr4(fiy_S, gmx_load_pr4(f+sciy)));
+        gmx_simd4_store_r(f+sciy, gmx_simd4_add_r(fiy_S, gmx_simd4_load_r(f+sciy)));
 
         fiz_S = gmx_mm_transpose_sum4h_pr(fiz_S0, fiz_S2);
-        gmx_store_pr4(f+sciz, gmx_add_pr4(fiz_S, gmx_load_pr4(f+sciz)));
+        gmx_simd4_store_r(f+sciz, gmx_simd4_add_r(fiz_S, gmx_simd4_load_r(f+sciz)));
 
 #ifdef CALC_SHIFTFORCES
-        fshift[ish3+0] += gmx_sum_simd4(fix_S, shf);
-        fshift[ish3+1] += gmx_sum_simd4(fiy_S, shf);
-        fshift[ish3+2] += gmx_sum_simd4(fiz_S, shf);
+        fshift[ish3+0] += gmx_simd4_reduce_r(fix_S);
+        fshift[ish3+1] += gmx_simd4_reduce_r(fiy_S);
+        fshift[ish3+2] += gmx_simd4_reduce_r(fiz_S);
 #endif
 
 #ifdef CALC_ENERGIES
         if (do_coul)
         {
-            *Vc += gmx_sum_simd(vctot_S, tmpsum);
+            *Vc += gmx_simd_reduce_r(vctot_S);
         }
-        *Vvdw += gmx_sum_simd(Vvdwtot_S, tmpsum);
+        *Vvdw += gmx_simd_reduce_r(Vvdwtot_S);
 #endif
 
         /* Outer loop uses 6 flops/iteration */