src/gromacs/nbnxm/kernels_simd_4xm/kernel_outer.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014,2015,2016 by the GROMACS development team.
   5  * Copyright (c) 2017,2018,2019,2020, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36
  37
  38 {
  39     using namespace gmx;
  40
  41     /* Unpack pointers for output */
  42     real* f      = out->f.data();
  43     real* fshift = out->fshift.data();
  44 #ifdef CALC_ENERGIES
  45 #    ifdef ENERGY_GROUPS
  46     real* Vvdw = out->VSvdw.data();
  47     real* Vc   = out->VSc.data();
  48 #    else
  49     real* Vvdw = out->Vvdw.data();
  50     real* Vc   = out->Vc.data();
  51 #    endif
  52 #endif
  53
  54     const nbnxn_cj_t* l_cj;
  55     int               ci, ci_sh;
  56     int               ish, ish3;
  57     gmx_bool          do_LJ, half_LJ, do_coul;
  58     int               cjind0, cjind1, cjind;
  59
  60 #ifdef ENERGY_GROUPS
  61     int   Vstride_i;
  62     int   egps_ishift, egps_imask;
  63     int   egps_jshift, egps_jmask, egps_jstride;
  64     int   egps_i;
  65     real* vvdwtp[UNROLLI];
  66     real* vctp[UNROLLI];
  67 #endif
  68
  69     SimdReal shX_S;
  70     SimdReal shY_S;
  71     SimdReal shZ_S;
  72     SimdReal ix_S0, iy_S0, iz_S0;
  73     SimdReal ix_S1, iy_S1, iz_S1;
  74     SimdReal ix_S2, iy_S2, iz_S2;
  75     SimdReal ix_S3, iy_S3, iz_S3;
  76     SimdReal fix_S0, fiy_S0, fiz_S0;
  77     SimdReal fix_S1, fiy_S1, fiz_S1;
  78     SimdReal fix_S2, fiy_S2, fiz_S2;
  79     SimdReal fix_S3, fiy_S3, fiz_S3;
  80
  81     SimdReal diagonal_jmi_S;
  82 #if UNROLLI == UNROLLJ
  83     SimdBool diagonal_mask_S0, diagonal_mask_S1, diagonal_mask_S2, diagonal_mask_S3;
  84 #else
  85     SimdBool diagonal_mask0_S0, diagonal_mask0_S1, diagonal_mask0_S2, diagonal_mask0_S3;
  86     SimdBool diagonal_mask1_S0, diagonal_mask1_S1, diagonal_mask1_S2, diagonal_mask1_S3;
  87 #endif
  88
  89     SimdBitMask filter_S0, filter_S1, filter_S2, filter_S3;
  90
  91     SimdReal zero_S(0.0);
  92
  93     SimdReal one_S(1.0);
  94     SimdReal iq_S0 = setZero();
  95     SimdReal iq_S1 = setZero();
  96     SimdReal iq_S2 = setZero();
  97     SimdReal iq_S3 = setZero();
  98
  99 #ifdef CALC_COUL_RF
 100     SimdReal mrc_3_S;
 101 #    ifdef CALC_ENERGIES
 102     SimdReal hrc_3_S, moh_rc_S;
 103 #    endif
 104 #endif
 105
 106 #ifdef CALC_COUL_TAB
 107     /* Coulomb table variables */
 108     SimdReal    invtsp_S;
 109     const real* tab_coul_F;
 110 #    if defined CALC_ENERGIES && !defined TAB_FDV0
 111     const real gmx_unused* tab_coul_V;
 112 #    endif
 113 #    ifdef CALC_ENERGIES
 114     SimdReal mhalfsp_S;
 115 #    endif
 116 #endif
 117
 118 #ifdef CALC_COUL_EWALD
 119     SimdReal beta2_S, beta_S;
 120 #endif
 121
 122 #if defined CALC_ENERGIES && (defined CALC_COUL_EWALD || defined CALC_COUL_TAB)
 123     SimdReal sh_ewald_S;
 124 #endif
 125
 126 #if defined LJ_CUT && defined CALC_ENERGIES
 127     SimdReal p6_cpot_S, p12_cpot_S;
 128 #endif
 129 #ifdef LJ_POT_SWITCH
 130     SimdReal rswitch_S;
 131     SimdReal swV3_S, swV4_S, swV5_S;
 132     SimdReal swF2_S, swF3_S, swF4_S;
 133 #endif
 134 #ifdef LJ_FORCE_SWITCH
 135     SimdReal rswitch_S;
 136     SimdReal p6_fc2_S, p6_fc3_S;
 137     SimdReal p12_fc2_S, p12_fc3_S;
 138 #    ifdef CALC_ENERGIES
 139     SimdReal p6_vc3_S, p6_vc4_S;
 140     SimdReal p12_vc3_S, p12_vc4_S;
 141     SimdReal p6_6cpot_S, p12_12cpot_S;
 142 #    endif
 143 #endif
 144 #ifdef LJ_EWALD_GEOM
 145     real     lj_ewaldcoeff2, lj_ewaldcoeff6_6;
 146     SimdReal mone_S, half_S, lje_c2_S, lje_c6_6_S;
 147 #endif
 148
 149 #ifdef LJ_COMB_LB
 150     SimdReal hsig_i_S0, seps_i_S0;
 151     SimdReal hsig_i_S1, seps_i_S1;
 152     SimdReal hsig_i_S2, seps_i_S2;
 153     SimdReal hsig_i_S3, seps_i_S3;
 154 #endif /* LJ_COMB_LB */
 155
 156     SimdReal minRsq_S;
 157     SimdReal rc2_S;
 158 #ifdef VDW_CUTOFF_CHECK
 159     SimdReal rcvdw2_S;
 160 #endif
 161
 162     int ninner;
 163
 164 #ifdef COUNT_PAIRS
 165     int npair = 0;
 166 #endif
 167
 168     const nbnxn_atomdata_t::Params& nbatParams = nbat->params();
 169
 170 #if defined LJ_COMB_GEOM || defined LJ_COMB_LB || defined LJ_EWALD_GEOM
 171     const real* gmx_restrict ljc = nbatParams.lj_comb.data();
 172 #endif
 173 #if !(defined LJ_COMB_GEOM || defined LJ_COMB_LB || defined FIX_LJ_C)
 174     /* No combination rule used */
 175     const real* gmx_restrict nbfp_ptr = nbatParams.nbfp_aligned.data();
 176     const int* gmx_restrict type      = nbatParams.type.data();
 177 #endif
 178
 179     /* Load j-i for the first i */
 180     diagonal_jmi_S = load<SimdReal>(nbat->simdMasks.diagonal_4xn_j_minus_i.data());
 181     /* Generate all the diagonal masks as comparison results */
 182 #if UNROLLI == UNROLLJ
 183     diagonal_mask_S0 = (zero_S < diagonal_jmi_S);
 184     diagonal_jmi_S   = diagonal_jmi_S - one_S;
 185     diagonal_mask_S1 = (zero_S < diagonal_jmi_S);
 186     diagonal_jmi_S   = diagonal_jmi_S - one_S;
 187     diagonal_mask_S2 = (zero_S < diagonal_jmi_S);
 188     diagonal_jmi_S   = diagonal_jmi_S - one_S;
 189     diagonal_mask_S3 = (zero_S < diagonal_jmi_S);
 190 #else
 191 #    if UNROLLI == 2 * UNROLLJ || 2 * UNROLLI == UNROLLJ
 192     diagonal_mask0_S0 = (zero_S < diagonal_jmi_S);
 193     diagonal_jmi_S    = diagonal_jmi_S - one_S;
 194     diagonal_mask0_S1 = (zero_S < diagonal_jmi_S);
 195     diagonal_jmi_S    = diagonal_jmi_S - one_S;
 196     diagonal_mask0_S2 = (zero_S < diagonal_jmi_S);
 197     diagonal_jmi_S    = diagonal_jmi_S - one_S;
 198     diagonal_mask0_S3 = (zero_S < diagonal_jmi_S);
 199     diagonal_jmi_S    = diagonal_jmi_S - one_S;
 200
 201 #        if UNROLLI == 2 * UNROLLJ
 202     /* Load j-i for the second half of the j-cluster */
 203     diagonal_jmi_S = load<SimdReal>(nbat->simdMasks.diagonal_4xn_j_minus_i.data() + UNROLLJ);
 204 #        endif
 205
 206     diagonal_mask1_S0                                  = (zero_S < diagonal_jmi_S);
 207     diagonal_jmi_S                                     = diagonal_jmi_S - one_S;
 208     diagonal_mask1_S1                                  = (zero_S < diagonal_jmi_S);
 209     diagonal_jmi_S                                     = diagonal_jmi_S - one_S;
 210     diagonal_mask1_S2                                  = (zero_S < diagonal_jmi_S);
 211     diagonal_jmi_S                                     = diagonal_jmi_S - one_S;
 212     diagonal_mask1_S3                                  = (zero_S < diagonal_jmi_S);
 213 #    endif
 214 #endif
 215
 216 #if GMX_DOUBLE && !GMX_SIMD_HAVE_INT32_LOGICAL
 217     const std::uint64_t* gmx_restrict exclusion_filter = nbat->simdMasks.exclusion_filter64.data();
 218 #else
 219     const std::uint32_t* gmx_restrict exclusion_filter = nbat->simdMasks.exclusion_filter.data();
 220 #endif
 221
 222     /* Here we cast the exclusion filters from unsigned * to int * or real *.
 223      * Since we only check bits, the actual value they represent does not
 224      * matter, as long as both filter and mask data are treated the same way.
 225      */
 226 #if GMX_SIMD_HAVE_INT32_LOGICAL
 227     filter_S0 = load<SimdBitMask>(reinterpret_cast<const int*>(exclusion_filter + 0 * UNROLLJ));
 228     filter_S1 = load<SimdBitMask>(reinterpret_cast<const int*>(exclusion_filter + 1 * UNROLLJ));
 229     filter_S2 = load<SimdBitMask>(reinterpret_cast<const int*>(exclusion_filter + 2 * UNROLLJ));
 230     filter_S3 = load<SimdBitMask>(reinterpret_cast<const int*>(exclusion_filter + 3 * UNROLLJ));
 231 #else
 232     filter_S0 = load<SimdBitMask>(reinterpret_cast<const real*>(exclusion_filter + 0 * UNROLLJ));
 233     filter_S1 = load<SimdBitMask>(reinterpret_cast<const real*>(exclusion_filter + 1 * UNROLLJ));
 234     filter_S2 = load<SimdBitMask>(reinterpret_cast<const real*>(exclusion_filter + 2 * UNROLLJ));
 235     filter_S3 = load<SimdBitMask>(reinterpret_cast<const real*>(exclusion_filter + 3 * UNROLLJ));
 236 #endif
 237
 238 #ifdef CALC_COUL_RF
 239     /* Reaction-field constants */
 240     mrc_3_S = SimdReal(-2 * ic->k_rf);
 241 #    ifdef CALC_ENERGIES
 242     hrc_3_S  = SimdReal(ic->k_rf);
 243     moh_rc_S = SimdReal(-ic->c_rf);
 244 #    endif
 245 #endif
 246
 247 #ifdef CALC_COUL_TAB
 248
 249     invtsp_S = SimdReal(ic->coulombEwaldTables->scale);
 250 #    ifdef CALC_ENERGIES
 251     mhalfsp_S = SimdReal(-0.5_real / ic->coulombEwaldTables->scale);
 252 #    endif
 253
 254 #    ifdef TAB_FDV0
 255     tab_coul_F = ic->coulombEwaldTables->tableFDV0.data();
 256 #    else
 257     tab_coul_F = ic->coulombEwaldTables->tableF.data();
 258 #        ifdef CALC_ENERGIES
 259     tab_coul_V = ic->coulombEwaldTables->tableV.data();
 260 #        endif
 261 #    endif
 262 #endif /* CALC_COUL_TAB */
 263
 264 #ifdef CALC_COUL_EWALD
 265     beta2_S = SimdReal(ic->ewaldcoeff_q * ic->ewaldcoeff_q);
 266     beta_S  = SimdReal(ic->ewaldcoeff_q);
 267 #endif
 268
 269 #if (defined CALC_COUL_TAB || defined CALC_COUL_EWALD) && defined CALC_ENERGIES
 270     sh_ewald_S = SimdReal(ic->sh_ewald);
 271 #endif
 272
 273     /* LJ function constants */
 274 #if defined CALC_ENERGIES || defined LJ_POT_SWITCH
 275     SimdReal sixth_S(1.0 / 6.0);
 276     SimdReal twelveth_S(1.0 / 12.0);
 277 #endif
 278
 279 #if defined LJ_CUT && defined CALC_ENERGIES
 280     /* We shift the potential by cpot, which can be zero */
 281     p6_cpot_S  = SimdReal(ic->dispersion_shift.cpot);
 282     p12_cpot_S = SimdReal(ic->repulsion_shift.cpot);
 283 #endif
 284 #ifdef LJ_POT_SWITCH
 285     rswitch_S = SimdReal(ic->rvdw_switch);
 286     swV3_S    = SimdReal(ic->vdw_switch.c3);
 287     swV4_S    = SimdReal(ic->vdw_switch.c4);
 288     swV5_S    = SimdReal(ic->vdw_switch.c5);
 289     swF2_S    = SimdReal(3 * ic->vdw_switch.c3);
 290     swF3_S    = SimdReal(4 * ic->vdw_switch.c4);
 291     swF4_S    = SimdReal(5 * ic->vdw_switch.c5);
 292 #endif
 293 #ifdef LJ_FORCE_SWITCH
 294     rswitch_S = SimdReal(ic->rvdw_switch);
 295     p6_fc2_S  = SimdReal(ic->dispersion_shift.c2);
 296     p6_fc3_S  = SimdReal(ic->dispersion_shift.c3);
 297     p12_fc2_S = SimdReal(ic->repulsion_shift.c2);
 298     p12_fc3_S = SimdReal(ic->repulsion_shift.c3);
 299 #    ifdef CALC_ENERGIES
 300     {
 301         SimdReal mthird_S(-1.0 / 3.0);
 302         SimdReal mfourth_S(-1.0 / 4.0);
 303
 304         p6_vc3_S     = mthird_S * p6_fc2_S;
 305         p6_vc4_S     = mfourth_S * p6_fc3_S;
 306         p6_6cpot_S   = SimdReal(ic->dispersion_shift.cpot / 6);
 307         p12_vc3_S    = mthird_S * p12_fc2_S;
 308         p12_vc4_S    = mfourth_S * p12_fc3_S;
 309         p12_12cpot_S = SimdReal(ic->repulsion_shift.cpot / 12);
 310     }
 311 #    endif
 312 #endif
 313 #ifdef LJ_EWALD_GEOM
 314     mone_S           = SimdReal(-1.0);
 315     half_S           = SimdReal(0.5);
 316     lj_ewaldcoeff2   = ic->ewaldcoeff_lj * ic->ewaldcoeff_lj;
 317     lj_ewaldcoeff6_6 = lj_ewaldcoeff2 * lj_ewaldcoeff2 * lj_ewaldcoeff2 / 6;
 318     lje_c2_S         = SimdReal(lj_ewaldcoeff2);
 319     lje_c6_6_S       = SimdReal(lj_ewaldcoeff6_6);
 320 #    ifdef CALC_ENERGIES
 321     /* Determine the grid potential at the cut-off */
 322     SimdReal lje_vc_S(ic->sh_lj_ewald);
 323 #    endif
 324 #endif
 325
 326     /* The kernel either supports rcoulomb = rvdw or rcoulomb >= rvdw */
 327     rc2_S = SimdReal(ic->rcoulomb * ic->rcoulomb);
 328 #ifdef VDW_CUTOFF_CHECK
 329     rcvdw2_S = SimdReal(ic->rvdw * ic->rvdw);
 330 #endif
 331
 332     minRsq_S = SimdReal(NBNXN_MIN_RSQ);
 333
 334     const real* gmx_restrict q        = nbatParams.q.data();
 335     const real               facel    = ic->epsfac;
 336     const real* gmx_restrict shiftvec = shift_vec[0];
 337     const real* gmx_restrict x        = nbat->x().data();
 338
 339 #ifdef FIX_LJ_C
 340     alignas(GMX_SIMD_ALIGNMENT) real pvdw_c6[2 * UNROLLI * UNROLLJ];
 341     real*                            pvdw_c12 = pvdw_c6 + UNROLLI * UNROLLJ;
 342
 343     for (int jp = 0; jp < UNROLLJ; jp++)
 344     {
 345         pvdw_c6[0 * UNROLLJ + jp] = nbat->nbfp[0 * 2];
 346         pvdw_c6[1 * UNROLLJ + jp] = nbat->nbfp[0 * 2];
 347         pvdw_c6[2 * UNROLLJ + jp] = nbat->nbfp[0 * 2];
 348         pvdw_c6[3 * UNROLLJ + jp] = nbat->nbfp[0 * 2];
 349
 350         pvdw_c12[0 * UNROLLJ + jp] = nbat->nbfp[0 * 2 + 1];
 351         pvdw_c12[1 * UNROLLJ + jp] = nbat->nbfp[0 * 2 + 1];
 352         pvdw_c12[2 * UNROLLJ + jp] = nbat->nbfp[0 * 2 + 1];
 353         pvdw_c12[3 * UNROLLJ + jp] = nbat->nbfp[0 * 2 + 1];
 354     }
 355     SimdReal c6_S0 = load<SimdReal>(pvdw_c6 + 0 * UNROLLJ);
 356     SimdReal c6_S1 = load<SimdReal>(pvdw_c6 + 1 * UNROLLJ);
 357     SimdReal c6_S2 = load<SimdReal>(pvdw_c6 + 2 * UNROLLJ);
 358     SimdReal c6_S3 = load<SimdReal>(pvdw_c6 + 3 * UNROLLJ);
 359
 360     SimdReal c12_S0 = load<SimdReal>(pvdw_c12 + 0 * UNROLLJ);
 361     SimdReal c12_S1 = load<SimdReal>(pvdw_c12 + 1 * UNROLLJ);
 362     SimdReal c12_S2 = load<SimdReal>(pvdw_c12 + 2 * UNROLLJ);
 363     SimdReal c12_S3 = load<SimdReal>(pvdw_c12 + 3 * UNROLLJ);
 364 #endif /* FIX_LJ_C */
 365
 366 #ifdef ENERGY_GROUPS
 367     egps_ishift  = nbatParams.neg_2log;
 368     egps_imask   = (1 << egps_ishift) - 1;
 369     egps_jshift  = 2 * nbatParams.neg_2log;
 370     egps_jmask   = (1 << egps_jshift) - 1;
 371     egps_jstride = (UNROLLJ >> 1) * UNROLLJ;
 372     /* Major division is over i-particle energy groups, determine the stride */
 373     Vstride_i = nbatParams.nenergrp * (1 << nbatParams.neg_2log) * egps_jstride;
 374 #endif
 375
 376     l_cj = nbl->cj.data();
 377
 378     ninner = 0;
 379
 380     for (const nbnxn_ci_t& ciEntry : nbl->ci)
 381     {
 382         ish    = (ciEntry.shift & NBNXN_CI_SHIFT);
 383         ish3   = ish * 3;
 384         cjind0 = ciEntry.cj_ind_start;
 385         cjind1 = ciEntry.cj_ind_end;
 386         ci     = ciEntry.ci;
 387         ci_sh  = (ish == CENTRAL ? ci : -1);
 388
 389         shX_S = SimdReal(shiftvec[ish3]);
 390         shY_S = SimdReal(shiftvec[ish3 + 1]);
 391         shZ_S = SimdReal(shiftvec[ish3 + 2]);
 392
 393 #if UNROLLJ <= 4
 394         int sci  = ci * STRIDE;
 395         int scix = sci * DIM;
 396 #    if defined LJ_COMB_LB || defined LJ_COMB_GEOM || defined LJ_EWALD_GEOM
 397         int sci2 = sci * 2;
 398 #    endif
 399 #else
 400         int sci  = (ci >> 1) * STRIDE;
 401         int scix = sci * DIM + (ci & 1) * (STRIDE >> 1);
 402 #    if defined LJ_COMB_LB || defined LJ_COMB_GEOM || defined LJ_EWALD_GEOM
 403         int sci2 = sci * 2 + (ci & 1) * (STRIDE >> 1);
 404 #    endif
 405         sci += (ci & 1) * (STRIDE >> 1);
 406 #endif
 407
 408         /* We have 5 LJ/C combinations, but use only three inner loops,
 409          * as the other combinations are unlikely and/or not much faster:
 410          * inner half-LJ + C for half-LJ + C / no-LJ + C
 411          * inner LJ + C      for full-LJ + C
 412          * inner LJ          for full-LJ + no-C / half-LJ + no-C
 413          */
 414         do_LJ   = ((ciEntry.shift & NBNXN_CI_DO_LJ(0)) != 0);
 415         do_coul = ((ciEntry.shift & NBNXN_CI_DO_COUL(0)) != 0);
 416         half_LJ = (((ciEntry.shift & NBNXN_CI_HALF_LJ(0)) != 0) || !do_LJ) && do_coul;
 417
 418 #ifdef ENERGY_GROUPS
 419         egps_i = nbatParams.energrp[ci];
 420         {
 421             int ia, egp_ia;
 422
 423             for (ia = 0; ia < UNROLLI; ia++)
 424             {
 425                 egp_ia     = (egps_i >> (ia * egps_ishift)) & egps_imask;
 426                 vvdwtp[ia] = Vvdw + egp_ia * Vstride_i;
 427                 vctp[ia]   = Vc + egp_ia * Vstride_i;
 428             }
 429         }
 430 #endif
 431
 432 #ifdef CALC_ENERGIES
 433 #    ifdef LJ_EWALD_GEOM
 434         gmx_bool do_self = TRUE;
 435 #    else
 436         gmx_bool do_self = do_coul;
 437 #    endif
 438 #    if UNROLLJ == 4
 439         if (do_self && l_cj[ciEntry.cj_ind_start].cj == ci_sh)
 440 #    endif
 441 #    if UNROLLJ == 2
 442             if (do_self && l_cj[ciEntry.cj_ind_start].cj == (ci_sh << 1))
 443 #    endif
 444 #    if UNROLLJ == 8
 445                 if (do_self && l_cj[ciEntry.cj_ind_start].cj == (ci_sh >> 1))
 446 #    endif
 447                 {
 448                     if (do_coul)
 449                     {
 450                         real Vc_sub_self;
 451                         int  ia;
 452
 453 #    ifdef CALC_COUL_RF
 454                         Vc_sub_self = 0.5 * ic->c_rf;
 455 #    endif
 456 #    ifdef CALC_COUL_TAB
 457 #        ifdef TAB_FDV0
 458                         Vc_sub_self = 0.5 * tab_coul_F[2];
 459 #        else
 460                         Vc_sub_self = 0.5 * tab_coul_V[0];
 461 #        endif
 462 #    endif
 463 #    ifdef CALC_COUL_EWALD
 464                         /* beta/sqrt(pi) */
 465                         Vc_sub_self = 0.5 * ic->ewaldcoeff_q * M_2_SQRTPI;
 466 #    endif
 467
 468                         for (ia = 0; ia < UNROLLI; ia++)
 469                         {
 470                             real qi;
 471
 472                             qi = q[sci + ia];
 473 #    ifdef ENERGY_GROUPS
 474                             vctp[ia][((egps_i >> (ia * egps_ishift)) & egps_imask) * egps_jstride]
 475 #    else
 476                     Vc[0]
 477 #    endif
 478                                     -= facel * qi * qi * Vc_sub_self;
 479                         }
 480                     }
 481
 482 #    ifdef LJ_EWALD_GEOM
 483                     {
 484                         int ia;
 485
 486                         for (ia = 0; ia < UNROLLI; ia++)
 487                         {
 488                             real c6_i;
 489
 490                             c6_i = nbatParams.nbfp[nbatParams.type[sci + ia] * (nbatParams.numTypes + 1) * 2]
 491                                    / 6;
 492 #        ifdef ENERGY_GROUPS
 493                             vvdwtp[ia][((egps_i >> (ia * egps_ishift)) & egps_imask) * egps_jstride]
 494 #        else
 495                             Vvdw[0]
 496 #        endif
 497                                     += 0.5 * c6_i * lj_ewaldcoeff6_6;
 498                         }
 499                     }
 500 #    endif /* LJ_EWALD_GEOM */
 501                 }
 502 #endif
 503
 504         /* Load i atom data */
 505         int sciy = scix + STRIDE;
 506         int sciz = sciy + STRIDE;
 507         ix_S0    = SimdReal(x[scix]) + shX_S;
 508         ix_S1    = SimdReal(x[scix + 1]) + shX_S;
 509         ix_S2    = SimdReal(x[scix + 2]) + shX_S;
 510         ix_S3    = SimdReal(x[scix + 3]) + shX_S;
 511         iy_S0    = SimdReal(x[sciy]) + shY_S;
 512         iy_S1    = SimdReal(x[sciy + 1]) + shY_S;
 513         iy_S2    = SimdReal(x[sciy + 2]) + shY_S;
 514         iy_S3    = SimdReal(x[sciy + 3]) + shY_S;
 515         iz_S0    = SimdReal(x[sciz]) + shZ_S;
 516         iz_S1    = SimdReal(x[sciz + 1]) + shZ_S;
 517         iz_S2    = SimdReal(x[sciz + 2]) + shZ_S;
 518         iz_S3    = SimdReal(x[sciz + 3]) + shZ_S;
 519
 520         if (do_coul)
 521         {
 522             iq_S0 = SimdReal(facel * q[sci]);
 523             iq_S1 = SimdReal(facel * q[sci + 1]);
 524             iq_S2 = SimdReal(facel * q[sci + 2]);
 525             iq_S3 = SimdReal(facel * q[sci + 3]);
 526         }
 527
 528 #ifdef LJ_COMB_LB
 529         hsig_i_S0 = SimdReal(ljc[sci2 + 0]);
 530         hsig_i_S1 = SimdReal(ljc[sci2 + 1]);
 531         hsig_i_S2 = SimdReal(ljc[sci2 + 2]);
 532         hsig_i_S3 = SimdReal(ljc[sci2 + 3]);
 533         seps_i_S0 = SimdReal(ljc[sci2 + STRIDE + 0]);
 534         seps_i_S1 = SimdReal(ljc[sci2 + STRIDE + 1]);
 535         seps_i_S2 = SimdReal(ljc[sci2 + STRIDE + 2]);
 536         seps_i_S3 = SimdReal(ljc[sci2 + STRIDE + 3]);
 537 #else
 538 #    ifdef LJ_COMB_GEOM
 539         SimdReal c6s_S0 = SimdReal(ljc[sci2 + 0]);
 540         SimdReal c6s_S1 = SimdReal(ljc[sci2 + 1]);
 541         SimdReal c6s_S2, c6s_S3;
 542         if (!half_LJ)
 543         {
 544             c6s_S2 = SimdReal(ljc[sci2 + 2]);
 545             c6s_S3 = SimdReal(ljc[sci2 + 3]);
 546         }
 547         else
 548         {
 549             c6s_S2 = setZero();
 550             c6s_S3 = setZero();
 551         }
 552         SimdReal c12s_S0 = SimdReal(ljc[sci2 + STRIDE + 0]);
 553         SimdReal c12s_S1 = SimdReal(ljc[sci2 + STRIDE + 1]);
 554         SimdReal c12s_S2, c12s_S3;
 555         if (!half_LJ)
 556         {
 557             c12s_S2 = SimdReal(ljc[sci2 + STRIDE + 2]);
 558             c12s_S3 = SimdReal(ljc[sci2 + STRIDE + 3]);
 559         }
 560         else
 561         {
 562             c12s_S2 = setZero();
 563             c12s_S3 = setZero();
 564         }
 565 #    else
 566         const int   numTypes = nbatParams.numTypes;
 567         const real* nbfp0    = nbfp_ptr + type[sci] * numTypes * c_simdBestPairAlignment;
 568         const real* nbfp1    = nbfp_ptr + type[sci + 1] * numTypes * c_simdBestPairAlignment;
 569         const real *nbfp2 = nullptr, *nbfp3 = nullptr;
 570         if (!half_LJ)
 571         {
 572             nbfp2 = nbfp_ptr + type[sci + 2] * numTypes * c_simdBestPairAlignment;
 573             nbfp3 = nbfp_ptr + type[sci + 3] * numTypes * c_simdBestPairAlignment;
 574         }
 575 #    endif
 576 #endif
 577 #ifdef LJ_EWALD_GEOM
 578         /* We need the geometrically combined C6 for the PME grid correction */
 579         SimdReal c6s_S0 = SimdReal(ljc[sci2 + 0]);
 580         SimdReal c6s_S1 = SimdReal(ljc[sci2 + 1]);
 581         SimdReal c6s_S2, c6s_S3;
 582         if (!half_LJ)
 583         {
 584             c6s_S2 = SimdReal(ljc[sci2 + 2]);
 585             c6s_S3 = SimdReal(ljc[sci2 + 3]);
 586         }
 587 #endif
 588
 589         /* Zero the potential energy for this list */
 590 #ifdef CALC_ENERGIES
 591         SimdReal Vvdwtot_S = setZero();
 592         SimdReal vctot_S   = setZero();
 593 #endif
 594
 595         /* Clear i atom forces */
 596         fix_S0 = setZero();
 597         fix_S1 = setZero();
 598         fix_S2 = setZero();
 599         fix_S3 = setZero();
 600         fiy_S0 = setZero();
 601         fiy_S1 = setZero();
 602         fiy_S2 = setZero();
 603         fiy_S3 = setZero();
 604         fiz_S0 = setZero();
 605         fiz_S1 = setZero();
 606         fiz_S2 = setZero();
 607         fiz_S3 = setZero();
 608
 609         cjind = cjind0;
 610
 611         /* Currently all kernels use (at least half) LJ */
 612 #define CALC_LJ
 613         if (half_LJ)
 614         {
 615             /* Coulomb: all i-atoms, LJ: first half i-atoms */
 616 #define CALC_COULOMB
 617 #define HALF_LJ
 618 #define CHECK_EXCLS
 619             while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
 620             {
 621 #include "kernel_inner.h"
 622                 cjind++;
 623             }
 624 #undef CHECK_EXCLS
 625             for (; (cjind < cjind1); cjind++)
 626             {
 627 #include "kernel_inner.h"
 628             }
 629 #undef HALF_LJ
 630 #undef CALC_COULOMB
 631         }
 632         else if (do_coul)
 633         {
 634             /* Coulomb: all i-atoms, LJ: all i-atoms */
 635 #define CALC_COULOMB
 636 #define CHECK_EXCLS
 637             while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
 638             {
 639 #include "kernel_inner.h"
 640                 cjind++;
 641             }
 642 #undef CHECK_EXCLS
 643             for (; (cjind < cjind1); cjind++)
 644             {
 645 #include "kernel_inner.h"
 646             }
 647 #undef CALC_COULOMB
 648         }
 649         else
 650         {
 651             /* Coulomb: none, LJ: all i-atoms */
 652 #define CHECK_EXCLS
 653             while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
 654             {
 655 #include "kernel_inner.h"
 656                 cjind++;
 657             }
 658 #undef CHECK_EXCLS
 659             for (; (cjind < cjind1); cjind++)
 660             {
 661 #include "kernel_inner.h"
 662             }
 663         }
 664 #undef CALC_LJ
 665         ninner += cjind1 - cjind0;
 666
 667         /* Add accumulated i-forces to the force array */
 668         real fShiftX = reduceIncr4ReturnSum(f + scix, fix_S0, fix_S1, fix_S2, fix_S3);
 669         real fShiftY = reduceIncr4ReturnSum(f + sciy, fiy_S0, fiy_S1, fiy_S2, fiy_S3);
 670         real fShiftZ = reduceIncr4ReturnSum(f + sciz, fiz_S0, fiz_S1, fiz_S2, fiz_S3);
 671
 672
 673 #ifdef CALC_SHIFTFORCES
 674         fshift[ish3 + 0] += fShiftX;
 675         fshift[ish3 + 1] += fShiftY;
 676         fshift[ish3 + 2] += fShiftZ;
 677 #endif
 678
 679 #ifdef CALC_ENERGIES
 680         if (do_coul)
 681         {
 682             *Vc += reduce(vctot_S);
 683         }
 684
 685         *Vvdw += reduce(Vvdwtot_S);
 686 #endif
 687
 688         /* Outer loop uses 6 flops/iteration */
 689     }
 690
 691 #ifdef COUNT_PAIRS
 692     printf("atom pairs %d\n", npair);
 693 #endif
 694 }