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;
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;
+ gmx_simd4_real_t fix_S, fiy_S, fiz_S;
gmx_simd_real_t diagonal_jmi_S;
#if UNROLLI == UNROLLJ
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;
-#else
+#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_vc3_S, p12_vc4_S;
gmx_simd_real_t p6_6cpot_S, p12_12cpot_S;
#endif
-#else
-#ifdef CALC_ENERGIES
- gmx_simd_real_t p6_cpot_S, p12_cpot_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
gmx_simd_real_t c6_S2, c12_S2;
#endif
-#ifdef LJ_COMB_GEOM
+#if defined LJ_COMB_GEOM || defined LJ_EWALD_GEOM
const real *ljc;
gmx_simd_real_t c6s_S0, c12s_S0;
- gmx_simd_real_t c6s_S1, c12s_S1;
gmx_simd_real_t c6s_S2 = gmx_simd_setzero_r();
gmx_simd_real_t c12s_S2 = gmx_simd_setzero_r();
- gmx_simd_real_t c6s_S3 = gmx_simd_setzero_r();
- gmx_simd_real_t c12s_S3 = gmx_simd_setzero_r();
#endif
#endif /* LJ_COMB_LB */
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
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);
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);
-#else
- sixth_S = gmx_simd_set1_r(1.0/6.0);
- twelveth_S = gmx_simd_set1_r(1.0/12.0);
+#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);
p12_12cpot_S = gmx_simd_set1_r(ic->repulsion_shift.cpot/12);
}
#endif
-#else
- /* Plain LJ cut-off, with potential shift cpot, which can be 0 */
-#ifdef CALC_ENERGIES
- 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_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
-#endif /* LJ_POT_SWITCH */
/* The kernel either supports rcoulomb = rvdw or rcoulomb >= rvdw */
rc2_S = gmx_simd_set1_r(ic->rcoulomb*ic->rcoulomb);
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
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 */
/* Add accumulated i-forces to the force array */
fix_S = gmx_mm_transpose_sum4h_pr(fix_S0, fix_S2);
- gmx_simd4_store_r(f+scix, gmx_add_pr4(fix_S, gmx_load_pr4(f+scix)));
+ gmx_simd4_store_r(f+scix, gmx_add_pr4(fix_S, gmx_simd4_load_r(f+scix)));
fiy_S = gmx_mm_transpose_sum4h_pr(fiy_S0, fiy_S2);
- gmx_simd4_store_r(f+sciy, gmx_add_pr4(fiy_S, gmx_load_pr4(f+sciy)));
+ gmx_simd4_store_r(f+sciy, gmx_add_pr4(fiy_S, gmx_simd4_load_r(f+sciy)));
fiz_S = gmx_mm_transpose_sum4h_pr(fiz_S0, fiz_S2);
- gmx_simd4_store_r(f+sciz, gmx_add_pr4(fiz_S, gmx_load_pr4(f+sciz)));
+ gmx_simd4_store_r(f+sciz, gmx_add_pr4(fiz_S, gmx_simd4_load_r(f+sciz)));
#ifdef CALC_SHIFTFORCES
fshift[ish3+0] += gmx_sum_simd4(fix_S, shf);
{
*Vc += gmx_sum_simd(vctot_S, tmpsum);
}
+
*Vvdw += gmx_sum_simd(Vvdwtot_S, tmpsum);
#endif