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