src/gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_outer.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013, 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 /* Half-width SIMD operations are required here.
  38  * As the 4xn kernels are the "standard" kernels and some special operations
  39  * are required only here, we define those in nbnxn_kernel_simd_utils_...
  40  *
  41  * Half-width SIMD real type:
  42  * gmx_mm_hpr
  43  *
  44  * Half-width SIMD operations
  45  * Load reals at half-width aligned pointer b into half-width SIMD register a:
  46  * gmx_load_hpr(a, b)
  47  * Set all entries in half-width SIMD register *a to b:
  48  * gmx_set1_hpr(a, b)
  49  * Load one real at b and one real at b+1 into halves of a, respectively:
  50  * gmx_load1p1_pr(a, b)
  51  * Load reals at half-width aligned pointer b into two halves of a:
  52  * gmx_loaddh_pr(a, b)
  53  * Store half-width SIMD register b into half width aligned memory a:
  54  * gmx_store_hpr(a, b)
  55  * gmx_add_hpr(a, b)
  56  * gmx_sub_hpr(a, b)
  57  * Sum over 4 half SIMD registers:
  58  * gmx_sum4_hpr(a, b)
  59  * Sum the elements of halfs of each input register and store sums in out:
  60  * gmx_mm_transpose_sum4h_pr(a, b)
  61  * Extract two half-width registers *b, *c from a full width register a:
  62  * gmx_pr_to_2hpr(a, b, c)
  63  */
  64
  65
  66 {
  67     const nbnxn_ci_t   *nbln;
  68     const nbnxn_cj_t   *l_cj;
  69     const int          *type;
  70     const real         *q;
  71     const real         *shiftvec;
  72     const real         *x;
  73     const real         *nbfp0, *nbfp1, *nbfp2 = NULL, *nbfp3 = NULL;
  74     real                facel;
  75     real               *nbfp_ptr;
  76     int                 n, ci, ci_sh;
  77     int                 ish, ish3;
  78     gmx_bool            do_LJ, half_LJ, do_coul;
  79     int                 sci, scix, sciy, sciz, sci2;
  80     int                 cjind0, cjind1, cjind;
  81     int                 ip, jp;
  82
  83 #ifdef ENERGY_GROUPS
  84     int         Vstride_i;
  85     int         egps_ishift, egps_imask;
  86     int         egps_jshift, egps_jmask, egps_jstride;
  87     int         egps_i;
  88     real       *vvdwtp[UNROLLI];
  89     real       *vctp[UNROLLI];
  90 #endif
  91
  92     gmx_mm_pr  shX_S;
  93     gmx_mm_pr  shY_S;
  94     gmx_mm_pr  shZ_S;
  95     gmx_mm_pr  ix_S0, iy_S0, iz_S0;
  96     gmx_mm_pr  ix_S2, iy_S2, iz_S2;
  97     gmx_mm_pr  fix_S0, fiy_S0, fiz_S0;
  98     gmx_mm_pr  fix_S2, fiy_S2, fiz_S2;
  99     /* We use an i-force SIMD register width of 4 */
 100     /* The pr4 stuff is defined in nbnxn_kernel_simd_utils.h */
 101     gmx_mm_pr4 fix_S, fiy_S, fiz_S;
 102
 103     gmx_mm_pr  diagonal_jmi_S;
 104 #if UNROLLI == UNROLLJ
 105     gmx_mm_pb  diagonal_mask_S0, diagonal_mask_S2;
 106 #else
 107     gmx_mm_pb  diagonal_mask0_S0, diagonal_mask0_S2;
 108     gmx_mm_pb  diagonal_mask1_S0, diagonal_mask1_S2;
 109 #endif
 110
 111     unsigned      *exclusion_filter;
 112     gmx_exclfilter filter_S0, filter_S2;
 113
 114     gmx_mm_pr      zero_S = gmx_set1_pr(0);
 115
 116     gmx_mm_pr      one_S = gmx_set1_pr(1.0);
 117     gmx_mm_pr      iq_S0 = gmx_setzero_pr();
 118     gmx_mm_pr      iq_S2 = gmx_setzero_pr();
 119     gmx_mm_pr      mrc_3_S;
 120 #ifdef CALC_ENERGIES
 121     gmx_mm_pr      hrc_3_S, moh_rc_S;
 122 #endif
 123
 124 #ifdef CALC_COUL_TAB
 125     /* Coulomb table variables */
 126     gmx_mm_pr   invtsp_S;
 127     const real *tab_coul_F;
 128 #ifndef TAB_FDV0
 129     const real *tab_coul_V;
 130 #endif
 131     int         ti0_array[2*GMX_SIMD_WIDTH_HERE], *ti0;
 132     int         ti2_array[2*GMX_SIMD_WIDTH_HERE], *ti2;
 133 #ifdef CALC_ENERGIES
 134     gmx_mm_pr   mhalfsp_S;
 135 #endif
 136 #endif
 137
 138 #ifdef CALC_COUL_EWALD
 139     gmx_mm_pr beta2_S, beta_S;
 140 #endif
 141
 142 #if defined CALC_ENERGIES && (defined CALC_COUL_EWALD || defined CALC_COUL_TAB)
 143     gmx_mm_pr  sh_ewald_S;
 144 #endif
 145
 146 #ifdef LJ_COMB_LB
 147     const real *ljc;
 148
 149     gmx_mm_pr   hsig_i_S0, seps_i_S0;
 150     gmx_mm_pr   hsig_i_S2, seps_i_S2;
 151 #else
 152 #ifdef FIX_LJ_C
 153     real        pvdw_array[2*UNROLLI*UNROLLJ+GMX_SIMD_WIDTH_HERE];
 154     real       *pvdw_c6, *pvdw_c12;
 155     gmx_mm_pr   c6_S0, c12_S0;
 156     gmx_mm_pr   c6_S2, c12_S2;
 157 #endif
 158
 159 #ifdef LJ_COMB_GEOM
 160     const real *ljc;
 161
 162     gmx_mm_pr   c6s_S0, c12s_S0;
 163     gmx_mm_pr   c6s_S1, c12s_S1;
 164     gmx_mm_pr   c6s_S2 = gmx_setzero_pr(), c12s_S2 = gmx_setzero_pr();
 165     gmx_mm_pr   c6s_S3 = gmx_setzero_pr(), c12s_S3 = gmx_setzero_pr();
 166 #endif
 167 #endif /* LJ_COMB_LB */
 168
 169     gmx_mm_pr  vctot_S, Vvdwtot_S;
 170     gmx_mm_pr  sixth_S, twelveth_S;
 171
 172     gmx_mm_pr  avoid_sing_S;
 173     gmx_mm_pr  rc2_S;
 174 #ifdef VDW_CUTOFF_CHECK
 175     gmx_mm_pr  rcvdw2_S;
 176 #endif
 177
 178 #ifdef CALC_ENERGIES
 179     gmx_mm_pr  sh_invrc6_S, sh_invrc12_S;
 180
 181     /* cppcheck-suppress unassignedVariable */
 182     real       tmpsum_array[2*GMX_SIMD_WIDTH_HERE], *tmpsum;
 183 #endif
 184 #ifdef CALC_SHIFTFORCES
 185     /* cppcheck-suppress unassignedVariable */
 186     real       shf_array[2*GMX_SIMD_WIDTH_HERE], *shf;
 187 #endif
 188
 189     int ninner;
 190
 191 #ifdef COUNT_PAIRS
 192     int npair = 0;
 193 #endif
 194
 195 #if defined LJ_COMB_GEOM || defined LJ_COMB_LB
 196     ljc = nbat->lj_comb;
 197 #else
 198     /* No combination rule used */
 199     nbfp_ptr    = (4 == nbfp_stride) ? nbat->nbfp_s4 : nbat->nbfp;
 200 #endif
 201
 202     /* Load j-i for the first i */
 203     diagonal_jmi_S    = gmx_load_pr(nbat->simd_2xnn_diagonal_j_minus_i);
 204     /* Generate all the diagonal masks as comparison results */
 205 #if UNROLLI == UNROLLJ
 206     diagonal_mask_S0  = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
 207     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
 208     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
 209     diagonal_mask_S2  = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
 210 #else
 211 #if 2*UNROLLI == UNROLLJ
 212     diagonal_mask0_S0 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
 213     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
 214     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
 215     diagonal_mask0_S2 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
 216     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
 217     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
 218     diagonal_mask1_S0 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
 219     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
 220     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
 221     diagonal_mask1_S2 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
 222 #endif
 223 #endif
 224
 225     /* Load masks for topology exclusion masking. filter_stride is
 226        static const, so the conditional will be optimized away. */
 227     if (1 == filter_stride)
 228     {
 229         exclusion_filter = nbat->simd_exclusion_filter1;
 230     }
 231     else /* (2 == filter_stride) */
 232     {
 233         exclusion_filter = nbat->simd_exclusion_filter2;
 234     }
 235
 236     /* Here we cast the exclusion masks from unsigned * to int * or
 237      * real *.  Since we only check bits, the actual value they
 238      * represent does not matter, as long as both mask and exclusion
 239      * info are treated the same way.
 240      */
 241     filter_S0 = gmx_load_exclusion_filter(exclusion_filter + 0*2*UNROLLJ*filter_stride);
 242     filter_S2 = gmx_load_exclusion_filter(exclusion_filter + 1*2*UNROLLJ*filter_stride);
 243
 244 #ifdef CALC_COUL_TAB
 245     /* Generate aligned table index pointers */
 246     ti0 = prepare_table_load_buffer(ti0_array);
 247     ti2 = prepare_table_load_buffer(ti2_array);
 248
 249     invtsp_S  = gmx_set1_pr(ic->tabq_scale);
 250 #ifdef CALC_ENERGIES
 251     mhalfsp_S = gmx_set1_pr(-0.5/ic->tabq_scale);
 252 #endif
 253
 254 #ifdef TAB_FDV0
 255     tab_coul_F = ic->tabq_coul_FDV0;
 256 #else
 257     tab_coul_F = ic->tabq_coul_F;
 258     tab_coul_V = ic->tabq_coul_V;
 259 #endif
 260 #endif /* CALC_COUL_TAB */
 261
 262 #ifdef CALC_COUL_EWALD
 263     beta2_S = gmx_set1_pr(ic->ewaldcoeff_q*ic->ewaldcoeff_q);
 264     beta_S  = gmx_set1_pr(ic->ewaldcoeff_q);
 265 #endif
 266
 267 #if (defined CALC_COUL_TAB || defined CALC_COUL_EWALD) && defined CALC_ENERGIES
 268     sh_ewald_S = gmx_set1_pr(ic->sh_ewald);
 269 #endif
 270
 271     q                   = nbat->q;
 272     type                = nbat->type;
 273     facel               = ic->epsfac;
 274     shiftvec            = shift_vec[0];
 275     x                   = nbat->x;
 276
 277     avoid_sing_S = gmx_set1_pr(NBNXN_AVOID_SING_R2_INC);
 278
 279     /* The kernel either supports rcoulomb = rvdw or rcoulomb >= rvdw */
 280     rc2_S    = gmx_set1_pr(ic->rcoulomb*ic->rcoulomb);
 281 #ifdef VDW_CUTOFF_CHECK
 282     rcvdw2_S = gmx_set1_pr(ic->rvdw*ic->rvdw);
 283 #endif
 284
 285 #ifdef CALC_ENERGIES
 286     sixth_S      = gmx_set1_pr(1.0/6.0);
 287     twelveth_S   = gmx_set1_pr(1.0/12.0);
 288
 289     sh_invrc6_S  = gmx_set1_pr(ic->sh_invrc6);
 290     sh_invrc12_S = gmx_set1_pr(ic->sh_invrc6*ic->sh_invrc6);
 291 #endif
 292
 293     mrc_3_S  = gmx_set1_pr(-2*ic->k_rf);
 294
 295 #ifdef CALC_ENERGIES
 296     hrc_3_S  = gmx_set1_pr(ic->k_rf);
 297
 298     moh_rc_S = gmx_set1_pr(-ic->c_rf);
 299 #endif
 300
 301 #ifdef CALC_ENERGIES
 302     tmpsum   = gmx_simd_align_real(tmpsum_array);
 303 #endif
 304 #ifdef CALC_SHIFTFORCES
 305     shf      = gmx_simd_align_real(shf_array);
 306 #endif
 307
 308 #ifdef FIX_LJ_C
 309     pvdw_c6  = gmx_simd_align_real(pvdw_array);
 310     pvdw_c12 = pvdw_c6 + UNROLLI*UNROLLJ;
 311
 312     for (jp = 0; jp < UNROLLJ; jp++)
 313     {
 314         pvdw_c6 [0*UNROLLJ+jp] = nbat->nbfp[0*2];
 315         pvdw_c6 [1*UNROLLJ+jp] = nbat->nbfp[0*2];
 316         pvdw_c6 [2*UNROLLJ+jp] = nbat->nbfp[0*2];
 317         pvdw_c6 [3*UNROLLJ+jp] = nbat->nbfp[0*2];
 318
 319         pvdw_c12[0*UNROLLJ+jp] = nbat->nbfp[0*2+1];
 320         pvdw_c12[1*UNROLLJ+jp] = nbat->nbfp[0*2+1];
 321         pvdw_c12[2*UNROLLJ+jp] = nbat->nbfp[0*2+1];
 322         pvdw_c12[3*UNROLLJ+jp] = nbat->nbfp[0*2+1];
 323     }
 324     c6_S0            = gmx_load_pr(pvdw_c6 +0*UNROLLJ);
 325     c6_S1            = gmx_load_pr(pvdw_c6 +1*UNROLLJ);
 326     c6_S2            = gmx_load_pr(pvdw_c6 +2*UNROLLJ);
 327     c6_S3            = gmx_load_pr(pvdw_c6 +3*UNROLLJ);
 328
 329     c12_S0           = gmx_load_pr(pvdw_c12+0*UNROLLJ);
 330     c12_S1           = gmx_load_pr(pvdw_c12+1*UNROLLJ);
 331     c12_S2           = gmx_load_pr(pvdw_c12+2*UNROLLJ);
 332     c12_S3           = gmx_load_pr(pvdw_c12+3*UNROLLJ);
 333 #endif /* FIX_LJ_C */
 334
 335 #ifdef ENERGY_GROUPS
 336     egps_ishift  = nbat->neg_2log;
 337     egps_imask   = (1<<egps_ishift) - 1;
 338     egps_jshift  = 2*nbat->neg_2log;
 339     egps_jmask   = (1<<egps_jshift) - 1;
 340     egps_jstride = (UNROLLJ>>1)*UNROLLJ;
 341     /* Major division is over i-particle energy groups, determine the stride */
 342     Vstride_i    = nbat->nenergrp*(1<<nbat->neg_2log)*egps_jstride;
 343 #endif
 344
 345     l_cj = nbl->cj;
 346
 347     ninner = 0;
 348     for (n = 0; n < nbl->nci; n++)
 349     {
 350         nbln = &nbl->ci[n];
 351
 352         ish              = (nbln->shift & NBNXN_CI_SHIFT);
 353         ish3             = ish*3;
 354         cjind0           = nbln->cj_ind_start;
 355         cjind1           = nbln->cj_ind_end;
 356         ci               = nbln->ci;
 357         ci_sh            = (ish == CENTRAL ? ci : -1);
 358
 359         shX_S = gmx_load1_pr(shiftvec+ish3);
 360         shY_S = gmx_load1_pr(shiftvec+ish3+1);
 361         shZ_S = gmx_load1_pr(shiftvec+ish3+2);
 362
 363 #if UNROLLJ <= 4
 364         sci              = ci*STRIDE;
 365         scix             = sci*DIM;
 366         sci2             = sci*2;
 367 #else
 368         sci              = (ci>>1)*STRIDE;
 369         scix             = sci*DIM + (ci & 1)*(STRIDE>>1);
 370         sci2             = sci*2 + (ci & 1)*(STRIDE>>1);
 371         sci             += (ci & 1)*(STRIDE>>1);
 372 #endif
 373
 374         /* We have 5 LJ/C combinations, but use only three inner loops,
 375          * as the other combinations are unlikely and/or not much faster:
 376          * inner half-LJ + C for half-LJ + C / no-LJ + C
 377          * inner LJ + C      for full-LJ + C
 378          * inner LJ          for full-LJ + no-C / half-LJ + no-C
 379          */
 380         do_LJ   = (nbln->shift & NBNXN_CI_DO_LJ(0));
 381         do_coul = (nbln->shift & NBNXN_CI_DO_COUL(0));
 382         half_LJ = ((nbln->shift & NBNXN_CI_HALF_LJ(0)) || !do_LJ) && do_coul;
 383
 384 #ifdef ENERGY_GROUPS
 385         egps_i = nbat->energrp[ci];
 386         {
 387             int ia, egp_ia;
 388
 389             for (ia = 0; ia < UNROLLI; ia++)
 390             {
 391                 egp_ia     = (egps_i >> (ia*egps_ishift)) & egps_imask;
 392                 vvdwtp[ia] = Vvdw + egp_ia*Vstride_i;
 393                 vctp[ia]   = Vc   + egp_ia*Vstride_i;
 394             }
 395         }
 396 #endif
 397 #if defined CALC_ENERGIES
 398 #if UNROLLJ == 4
 399         if (do_coul && l_cj[nbln->cj_ind_start].cj == ci_sh)
 400 #endif
 401 #if UNROLLJ == 8
 402         if (do_coul && l_cj[nbln->cj_ind_start].cj == (ci_sh>>1))
 403 #endif
 404         {
 405             int  ia;
 406             real Vc_sub_self;
 407
 408 #ifdef CALC_COUL_RF
 409             Vc_sub_self = 0.5*ic->c_rf;
 410 #endif
 411 #ifdef CALC_COUL_TAB
 412 #ifdef TAB_FDV0
 413             Vc_sub_self = 0.5*tab_coul_F[2];
 414 #else
 415             Vc_sub_self = 0.5*tab_coul_V[0];
 416 #endif
 417 #endif
 418 #ifdef CALC_COUL_EWALD
 419             /* beta/sqrt(pi) */
 420             Vc_sub_self = 0.5*ic->ewaldcoeff_q*M_2_SQRTPI;
 421 #endif
 422
 423             for (ia = 0; ia < UNROLLI; ia++)
 424             {
 425                 real qi;
 426
 427                 qi = q[sci+ia];
 428 #ifdef ENERGY_GROUPS
 429                 vctp[ia][((egps_i>>(ia*egps_ishift)) & egps_imask)*egps_jstride]
 430 #else
 431                 Vc[0]
 432 #endif
 433                     -= facel*qi*qi*Vc_sub_self;
 434             }
 435         }
 436 #endif
 437
 438         /* Load i atom data */
 439         sciy             = scix + STRIDE;
 440         sciz             = sciy + STRIDE;
 441         gmx_load1p1_pr(&ix_S0, x+scix);
 442         gmx_load1p1_pr(&ix_S2, x+scix+2);
 443         gmx_load1p1_pr(&iy_S0, x+sciy);
 444         gmx_load1p1_pr(&iy_S2, x+sciy+2);
 445         gmx_load1p1_pr(&iz_S0, x+sciz);
 446         gmx_load1p1_pr(&iz_S2, x+sciz+2);
 447         ix_S0          = gmx_add_pr(ix_S0, shX_S);
 448         ix_S2          = gmx_add_pr(ix_S2, shX_S);
 449         iy_S0          = gmx_add_pr(iy_S0, shY_S);
 450         iy_S2          = gmx_add_pr(iy_S2, shY_S);
 451         iz_S0          = gmx_add_pr(iz_S0, shZ_S);
 452         iz_S2          = gmx_add_pr(iz_S2, shZ_S);
 453
 454         if (do_coul)
 455         {
 456             gmx_mm_pr facel_S;
 457
 458             facel_S    = gmx_set1_pr(facel);
 459
 460             gmx_load1p1_pr(&iq_S0, q+sci);
 461             gmx_load1p1_pr(&iq_S2, q+sci+2);
 462             iq_S0      = gmx_mul_pr(facel_S, iq_S0);
 463             iq_S2      = gmx_mul_pr(facel_S, iq_S2);
 464         }
 465
 466 #ifdef LJ_COMB_LB
 467         gmx_load1p1_pr(&hsig_i_S0, ljc+sci2+0);
 468         gmx_load1p1_pr(&hsig_i_S2, ljc+sci2+2);
 469         gmx_load1p1_pr(&seps_i_S0, ljc+sci2+STRIDE+0);
 470         gmx_load1p1_pr(&seps_i_S2, ljc+sci2+STRIDE+2);
 471 #else
 472 #ifdef LJ_COMB_GEOM
 473         gmx_load1p1_pr(&c6s_S0, ljc+sci2+0);
 474         if (!half_LJ)
 475         {
 476             gmx_load1p1_pr(&c6s_S2, ljc+sci2+2);
 477         }
 478         gmx_load1p1_pr(&c12s_S0, ljc+sci2+STRIDE+0);
 479         if (!half_LJ)
 480         {
 481             gmx_load1p1_pr(&c12s_S2, ljc+sci2+STRIDE+2);
 482         }
 483 #else
 484         nbfp0     = nbfp_ptr + type[sci  ]*nbat->ntype*nbfp_stride;
 485         nbfp1     = nbfp_ptr + type[sci+1]*nbat->ntype*nbfp_stride;
 486         if (!half_LJ)
 487         {
 488             nbfp2 = nbfp_ptr + type[sci+2]*nbat->ntype*nbfp_stride;
 489             nbfp3 = nbfp_ptr + type[sci+3]*nbat->ntype*nbfp_stride;
 490         }
 491 #endif
 492 #endif
 493
 494         /* Zero the potential energy for this list */
 495         Vvdwtot_S        = gmx_setzero_pr();
 496         vctot_S          = gmx_setzero_pr();
 497
 498         /* Clear i atom forces */
 499         fix_S0           = gmx_setzero_pr();
 500         fix_S2           = gmx_setzero_pr();
 501         fiy_S0           = gmx_setzero_pr();
 502         fiy_S2           = gmx_setzero_pr();
 503         fiz_S0           = gmx_setzero_pr();
 504         fiz_S2           = gmx_setzero_pr();
 505
 506         cjind = cjind0;
 507
 508         /* Currently all kernels use (at least half) LJ */
 509 #define CALC_LJ
 510         if (half_LJ)
 511         {
 512 #define CALC_COULOMB
 513 #define HALF_LJ
 514 #define CHECK_EXCLS
 515             while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
 516             {
 517 #include "nbnxn_kernel_simd_2xnn_inner.h"
 518                 cjind++;
 519             }
 520 #undef CHECK_EXCLS
 521             for (; (cjind < cjind1); cjind++)
 522             {
 523 #include "nbnxn_kernel_simd_2xnn_inner.h"
 524             }
 525 #undef HALF_LJ
 526 #undef CALC_COULOMB
 527         }
 528         else if (do_coul)
 529         {
 530 #define CALC_COULOMB
 531 #define CHECK_EXCLS
 532             while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
 533             {
 534 #include "nbnxn_kernel_simd_2xnn_inner.h"
 535                 cjind++;
 536             }
 537 #undef CHECK_EXCLS
 538             for (; (cjind < cjind1); cjind++)
 539             {
 540 #include "nbnxn_kernel_simd_2xnn_inner.h"
 541             }
 542 #undef CALC_COULOMB
 543         }
 544         else
 545         {
 546 #define CHECK_EXCLS
 547             while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
 548             {
 549 #include "nbnxn_kernel_simd_2xnn_inner.h"
 550                 cjind++;
 551             }
 552 #undef CHECK_EXCLS
 553             for (; (cjind < cjind1); cjind++)
 554             {
 555 #include "nbnxn_kernel_simd_2xnn_inner.h"
 556             }
 557         }
 558 #undef CALC_LJ
 559         ninner += cjind1 - cjind0;
 560
 561         /* Add accumulated i-forces to the force array */
 562         fix_S = gmx_mm_transpose_sum4h_pr(fix_S0, fix_S2);
 563         gmx_store_pr4(f+scix, gmx_add_pr4(fix_S, gmx_load_pr4(f+scix)));
 564
 565         fiy_S = gmx_mm_transpose_sum4h_pr(fiy_S0, fiy_S2);
 566         gmx_store_pr4(f+sciy, gmx_add_pr4(fiy_S, gmx_load_pr4(f+sciy)));
 567
 568         fiz_S = gmx_mm_transpose_sum4h_pr(fiz_S0, fiz_S2);
 569         gmx_store_pr4(f+sciz, gmx_add_pr4(fiz_S, gmx_load_pr4(f+sciz)));
 570
 571 #ifdef CALC_SHIFTFORCES
 572         fshift[ish3+0] += gmx_sum_simd4(fix_S, shf);
 573         fshift[ish3+1] += gmx_sum_simd4(fiy_S, shf);
 574         fshift[ish3+2] += gmx_sum_simd4(fiz_S, shf);
 575 #endif
 576
 577 #ifdef CALC_ENERGIES
 578         if (do_coul)
 579         {
 580             *Vc += gmx_sum_simd(vctot_S, tmpsum);
 581         }
 582         *Vvdw += gmx_sum_simd(Vvdwtot_S, tmpsum);
 583 #endif
 584
 585         /* Outer loop uses 6 flops/iteration */
 586     }
 587
 588 #ifdef COUNT_PAIRS
 589     printf("atom pairs %d\n", npair);
 590 #endif
 591 }