src/gromacs/mdlib/nbnxn_kernels/simd_4xn/nbnxn_kernel_simd_4xn_outer.h

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