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