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