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