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) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2009, The GROMACS Development Team
   6  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   7  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   8  * others, as listed in the AUTHORS file in the top-level source
   9  * directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37
  38 {
  39     const nbnxn_ci_t   *nbln;
  40     const nbnxn_cj_t   *l_cj;
  41     const int          *type;
  42     const real         *q;
  43     const real         *shiftvec;
  44     const real         *x;
  45     const real         *nbfp0, *nbfp1, *nbfp2 = NULL, *nbfp3 = NULL;
  46     real                facel;
  47     real               *nbfp_ptr;
  48     int                 n, ci, ci_sh;
  49     int                 ish, ish3;
  50     gmx_bool            do_LJ, half_LJ, do_coul;
  51     int                 sci, scix, sciy, sciz, sci2;
  52     int                 cjind0, cjind1, cjind;
  53     int                 ip, jp;
  54
  55 #ifdef ENERGY_GROUPS
  56     int         Vstride_i;
  57     int         egps_ishift, egps_imask;
  58     int         egps_jshift, egps_jmask, egps_jstride;
  59     int         egps_i;
  60     real       *vvdwtp[UNROLLI];
  61     real       *vctp[UNROLLI];
  62 #endif
  63
  64     gmx_mm_pr  shX_S;
  65     gmx_mm_pr  shY_S;
  66     gmx_mm_pr  shZ_S;
  67     gmx_mm_pr  ix_S0, iy_S0, iz_S0;
  68     gmx_mm_pr  ix_S1, iy_S1, iz_S1;
  69     gmx_mm_pr  ix_S2, iy_S2, iz_S2;
  70     gmx_mm_pr  ix_S3, iy_S3, iz_S3;
  71     gmx_mm_pr  fix_S0, fiy_S0, fiz_S0;
  72     gmx_mm_pr  fix_S1, fiy_S1, fiz_S1;
  73     gmx_mm_pr  fix_S2, fiy_S2, fiz_S2;
  74     gmx_mm_pr  fix_S3, fiy_S3, fiz_S3;
  75 #if UNROLLJ >= 4
  76     /* We use an i-force SIMD register width of 4 */
  77 #if UNROLLJ == 4
  78 #define gmx_mm_pr4     gmx_mm_pr
  79 #define gmx_load_pr4   gmx_load_pr
  80 #define gmx_store_pr4  gmx_store_pr
  81 #define gmx_add_pr4    gmx_add_pr
  82 #else
  83     /* The pr4 stuff is defined in nbnxn_kernel_simd_utils.h */
  84 #endif
  85     gmx_mm_pr4 fix_S, fiy_S, fiz_S;
  86 #else
  87     /* We use an i-force SIMD register width of 2 */
  88     gmx_mm_pr  fix0_S, fiy0_S, fiz0_S;
  89     gmx_mm_pr  fix2_S, fiy2_S, fiz2_S;
  90 #endif
  91
  92     gmx_mm_pr  diagonal_jmi_S;
  93 #if UNROLLI == UNROLLJ
  94     gmx_mm_pb  diagonal_mask_S0, diagonal_mask_S1, diagonal_mask_S2, diagonal_mask_S3;
  95 #else
  96     gmx_mm_pb  diagonal_mask0_S0, diagonal_mask0_S1, diagonal_mask0_S2, diagonal_mask0_S3;
  97     gmx_mm_pb  diagonal_mask1_S0, diagonal_mask1_S1, diagonal_mask1_S2, diagonal_mask1_S3;
  98 #endif
  99
 100     unsigned      *exclusion_filter;
 101     gmx_exclfilter filter_S0, filter_S1, filter_S2, filter_S3;
 102
 103     gmx_mm_pr  zero_S = gmx_set1_pr(0);
 104
 105     gmx_mm_pr  one_S  = gmx_set1_pr(1.0);
 106     gmx_mm_pr  iq_S0  = gmx_setzero_pr();
 107     gmx_mm_pr  iq_S1  = gmx_setzero_pr();
 108     gmx_mm_pr  iq_S2  = gmx_setzero_pr();
 109     gmx_mm_pr  iq_S3  = gmx_setzero_pr();
 110     gmx_mm_pr  mrc_3_S;
 111 #ifdef CALC_ENERGIES
 112     gmx_mm_pr  hrc_3_S, moh_rc_S;
 113 #endif
 114
 115 #ifdef CALC_COUL_TAB
 116     /* Coulomb table variables */
 117     gmx_mm_pr   invtsp_S;
 118     const real *tab_coul_F;
 119 #ifndef TAB_FDV0
 120     const real *tab_coul_V;
 121 #endif
 122     /* Thread-local working buffers for force and potential lookups */
 123     int        ti0_array[2*GMX_SIMD_WIDTH_HERE-1], *ti0 = NULL;
 124     int        ti1_array[2*GMX_SIMD_WIDTH_HERE-1], *ti1 = NULL;
 125     int        ti2_array[2*GMX_SIMD_WIDTH_HERE-1], *ti2 = NULL;
 126     int        ti3_array[2*GMX_SIMD_WIDTH_HERE-1], *ti3 = NULL;
 127 #ifdef CALC_ENERGIES
 128     gmx_mm_pr  mhalfsp_S;
 129 #endif
 130 #endif
 131
 132 #ifdef CALC_COUL_EWALD
 133     gmx_mm_pr beta2_S, beta_S;
 134 #endif
 135
 136 #if defined CALC_ENERGIES && (defined CALC_COUL_EWALD || defined CALC_COUL_TAB)
 137     gmx_mm_pr  sh_ewald_S;
 138 #endif
 139
 140 #ifdef LJ_COMB_LB
 141     const real *ljc;
 142
 143     gmx_mm_pr   hsig_i_S0, seps_i_S0;
 144     gmx_mm_pr   hsig_i_S1, seps_i_S1;
 145     gmx_mm_pr   hsig_i_S2, seps_i_S2;
 146     gmx_mm_pr   hsig_i_S3, seps_i_S3;
 147 #else
 148 #ifdef FIX_LJ_C
 149     real        pvdw_array[2*UNROLLI*UNROLLJ+3];
 150     real       *pvdw_c6, *pvdw_c12;
 151     gmx_mm_pr   c6_S0, c12_S0;
 152     gmx_mm_pr   c6_S1, c12_S1;
 153     gmx_mm_pr   c6_S2, c12_S2;
 154     gmx_mm_pr   c6_S3, c12_S3;
 155 #endif
 156
 157 #ifdef LJ_COMB_GEOM
 158     const real *ljc;
 159
 160     gmx_mm_pr   c6s_S0, c12s_S0;
 161     gmx_mm_pr   c6s_S1, c12s_S1;
 162     gmx_mm_pr   c6s_S2 = gmx_setzero_pr(), c12s_S2 = gmx_setzero_pr();
 163     gmx_mm_pr   c6s_S3 = gmx_setzero_pr(), c12s_S3 = gmx_setzero_pr();
 164 #endif
 165 #endif /* LJ_COMB_LB */
 166
 167     gmx_mm_pr  vctot_S, Vvdwtot_S;
 168     gmx_mm_pr  sixth_S, twelveth_S;
 169
 170     gmx_mm_pr  avoid_sing_S;
 171     gmx_mm_pr  rc2_S;
 172 #ifdef VDW_CUTOFF_CHECK
 173     gmx_mm_pr  rcvdw2_S;
 174 #endif
 175
 176 #ifdef CALC_ENERGIES
 177     gmx_mm_pr  sh_invrc6_S, sh_invrc12_S;
 178
 179     /* cppcheck-suppress unassignedVariable */
 180     real       tmpsum_array[15], *tmpsum;
 181 #endif
 182 #ifdef CALC_SHIFTFORCES
 183     /* cppcheck-suppress unassignedVariable */
 184     real       shf_array[15], *shf;
 185 #endif
 186
 187     int ninner;
 188
 189 #ifdef COUNT_PAIRS
 190     int npair = 0;
 191 #endif
 192
 193 #if defined LJ_COMB_GEOM || defined LJ_COMB_LB
 194     ljc = nbat->lj_comb;
 195 #else
 196     /* No combination rule used */
 197     nbfp_ptr    = (4 == nbfp_stride) ? nbat->nbfp_s4 : nbat->nbfp;
 198 #endif
 199
 200     /* Load j-i for the first i */
 201     diagonal_jmi_S    = gmx_load_pr(nbat->simd_4xn_diagonal_j_minus_i);
 202     /* Generate all the diagonal masks as comparison results */
 203 #if UNROLLI == UNROLLJ
 204     diagonal_mask_S0  = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
 205     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
 206     diagonal_mask_S1  = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
 207     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
 208     diagonal_mask_S2  = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
 209     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
 210     diagonal_mask_S3  = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
 211 #else
 212 #if UNROLLI == 2*UNROLLJ || 2*UNROLLI == UNROLLJ
 213     diagonal_mask0_S0 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
 214     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
 215     diagonal_mask0_S1 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
 216     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
 217     diagonal_mask0_S2 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
 218     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
 219     diagonal_mask0_S3 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
 220     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
 221
 222 #if UNROLLI == 2*UNROLLJ
 223     /* Load j-i for the second half of the j-cluster */
 224     diagonal_jmi_S    = gmx_load_pr(nbat->simd_4xn_diagonal_j_minus_i + UNROLLJ);
 225 #endif
 226
 227     diagonal_mask1_S0 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
 228     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
 229     diagonal_mask1_S1 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
 230     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
 231     diagonal_mask1_S2 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
 232     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
 233     diagonal_mask1_S3 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
 234 #endif
 235 #endif
 236
 237     /* Load masks for topology exclusion masking. filter_stride is
 238        static const, so the conditional will be optimized away. */
 239     if (1 == filter_stride)
 240     {
 241         exclusion_filter = nbat->simd_exclusion_filter1;
 242     }
 243     else /* (2 == filter_stride) */
 244     {
 245         exclusion_filter = nbat->simd_exclusion_filter2;
 246     }
 247
 248     /* Here we cast the exclusion filters from unsigned * to int * or real *.
 249      * Since we only check bits, the actual value they represent does not
 250      * matter, as long as both filter and mask data are treated the same way.
 251      */
 252     filter_S0    = gmx_load_exclusion_filter(exclusion_filter + 0*UNROLLJ*filter_stride);
 253     filter_S1    = gmx_load_exclusion_filter(exclusion_filter + 1*UNROLLJ*filter_stride);
 254     filter_S2    = gmx_load_exclusion_filter(exclusion_filter + 2*UNROLLJ*filter_stride);
 255     filter_S3    = gmx_load_exclusion_filter(exclusion_filter + 3*UNROLLJ*filter_stride);
 256
 257 #ifdef CALC_COUL_TAB
 258     /* Generate aligned table index pointers */
 259     ti0 = prepare_table_load_buffer(ti0_array);
 260     ti1 = prepare_table_load_buffer(ti1_array);
 261     ti2 = prepare_table_load_buffer(ti2_array);
 262     ti3 = prepare_table_load_buffer(ti3_array);
 263
 264     invtsp_S  = gmx_set1_pr(ic->tabq_scale);
 265 #ifdef CALC_ENERGIES
 266     mhalfsp_S = gmx_set1_pr(-0.5/ic->tabq_scale);
 267 #endif
 268
 269 #ifdef TAB_FDV0
 270     tab_coul_F = ic->tabq_coul_FDV0;
 271 #else
 272     tab_coul_F = ic->tabq_coul_F;
 273     tab_coul_V = ic->tabq_coul_V;
 274 #endif
 275 #endif /* CALC_COUL_TAB */
 276
 277 #ifdef CALC_COUL_EWALD
 278     beta2_S = gmx_set1_pr(ic->ewaldcoeff*ic->ewaldcoeff);
 279     beta_S  = gmx_set1_pr(ic->ewaldcoeff);
 280 #endif
 281
 282 #if (defined CALC_COUL_TAB || defined CALC_COUL_EWALD) && defined CALC_ENERGIES
 283     sh_ewald_S = gmx_set1_pr(ic->sh_ewald);
 284 #endif
 285
 286     q                   = nbat->q;
 287     type                = nbat->type;
 288     facel               = ic->epsfac;
 289     shiftvec            = shift_vec[0];
 290     x                   = nbat->x;
 291
 292     avoid_sing_S = gmx_set1_pr(NBNXN_AVOID_SING_R2_INC);
 293
 294     /* The kernel either supports rcoulomb = rvdw or rcoulomb >= rvdw */
 295     rc2_S    = gmx_set1_pr(ic->rcoulomb*ic->rcoulomb);
 296 #ifdef VDW_CUTOFF_CHECK
 297     rcvdw2_S = gmx_set1_pr(ic->rvdw*ic->rvdw);
 298 #endif
 299
 300 #ifdef CALC_ENERGIES
 301     sixth_S      = gmx_set1_pr(1.0/6.0);
 302     twelveth_S   = gmx_set1_pr(1.0/12.0);
 303
 304     sh_invrc6_S  = gmx_set1_pr(ic->sh_invrc6);
 305     sh_invrc12_S = gmx_set1_pr(ic->sh_invrc6*ic->sh_invrc6);
 306 #endif
 307
 308     mrc_3_S  = gmx_set1_pr(-2*ic->k_rf);
 309
 310 #ifdef CALC_ENERGIES
 311     hrc_3_S  = gmx_set1_pr(ic->k_rf);
 312
 313     moh_rc_S = gmx_set1_pr(-ic->c_rf);
 314 #endif
 315
 316 #ifdef CALC_ENERGIES
 317     tmpsum   = gmx_simd_align_real(tmpsum_array);
 318 #endif
 319 #ifdef CALC_SHIFTFORCES
 320     shf      = gmx_simd_align_real(shf_array);
 321 #endif
 322
 323 #ifdef FIX_LJ_C
 324     pvdw_c6  = gmx_simd_align_real(pvdw_array+3);
 325     pvdw_c12 = pvdw_c6 + UNROLLI*UNROLLJ;
 326
 327     for (jp = 0; jp < UNROLLJ; jp++)
 328     {
 329         pvdw_c6 [0*UNROLLJ+jp] = nbat->nbfp[0*2];
 330         pvdw_c6 [1*UNROLLJ+jp] = nbat->nbfp[0*2];
 331         pvdw_c6 [2*UNROLLJ+jp] = nbat->nbfp[0*2];
 332         pvdw_c6 [3*UNROLLJ+jp] = nbat->nbfp[0*2];
 333
 334         pvdw_c12[0*UNROLLJ+jp] = nbat->nbfp[0*2+1];
 335         pvdw_c12[1*UNROLLJ+jp] = nbat->nbfp[0*2+1];
 336         pvdw_c12[2*UNROLLJ+jp] = nbat->nbfp[0*2+1];
 337         pvdw_c12[3*UNROLLJ+jp] = nbat->nbfp[0*2+1];
 338     }
 339     c6_S0            = gmx_load_pr(pvdw_c6 +0*UNROLLJ);
 340     c6_S1            = gmx_load_pr(pvdw_c6 +1*UNROLLJ);
 341     c6_S2            = gmx_load_pr(pvdw_c6 +2*UNROLLJ);
 342     c6_S3            = gmx_load_pr(pvdw_c6 +3*UNROLLJ);
 343
 344     c12_S0           = gmx_load_pr(pvdw_c12+0*UNROLLJ);
 345     c12_S1           = gmx_load_pr(pvdw_c12+1*UNROLLJ);
 346     c12_S2           = gmx_load_pr(pvdw_c12+2*UNROLLJ);
 347     c12_S3           = gmx_load_pr(pvdw_c12+3*UNROLLJ);
 348 #endif /* FIX_LJ_C */
 349
 350 #ifdef ENERGY_GROUPS
 351     egps_ishift  = nbat->neg_2log;
 352     egps_imask   = (1<<egps_ishift) - 1;
 353     egps_jshift  = 2*nbat->neg_2log;
 354     egps_jmask   = (1<<egps_jshift) - 1;
 355     egps_jstride = (UNROLLJ>>1)*UNROLLJ;
 356     /* Major division is over i-particle energy groups, determine the stride */
 357     Vstride_i    = nbat->nenergrp*(1<<nbat->neg_2log)*egps_jstride;
 358 #endif
 359
 360     l_cj = nbl->cj;
 361
 362     ninner = 0;
 363     for (n = 0; n < nbl->nci; n++)
 364     {
 365         nbln = &nbl->ci[n];
 366
 367         ish              = (nbln->shift & NBNXN_CI_SHIFT);
 368         ish3             = ish*3;
 369         cjind0           = nbln->cj_ind_start;
 370         cjind1           = nbln->cj_ind_end;
 371         ci               = nbln->ci;
 372         ci_sh            = (ish == CENTRAL ? ci : -1);
 373
 374         shX_S = gmx_load1_pr(shiftvec+ish3);
 375         shY_S = gmx_load1_pr(shiftvec+ish3+1);
 376         shZ_S = gmx_load1_pr(shiftvec+ish3+2);
 377
 378 #if UNROLLJ <= 4
 379         sci              = ci*STRIDE;
 380         scix             = sci*DIM;
 381         sci2             = sci*2;
 382 #else
 383         sci              = (ci>>1)*STRIDE;
 384         scix             = sci*DIM + (ci & 1)*(STRIDE>>1);
 385         sci2             = sci*2 + (ci & 1)*(STRIDE>>1);
 386         sci             += (ci & 1)*(STRIDE>>1);
 387 #endif
 388
 389         /* We have 5 LJ/C combinations, but use only three inner loops,
 390          * as the other combinations are unlikely and/or not much faster:
 391          * inner half-LJ + C for half-LJ + C / no-LJ + C
 392          * inner LJ + C      for full-LJ + C
 393          * inner LJ          for full-LJ + no-C / half-LJ + no-C
 394          */
 395         do_LJ   = (nbln->shift & NBNXN_CI_DO_LJ(0));
 396         do_coul = (nbln->shift & NBNXN_CI_DO_COUL(0));
 397         half_LJ = ((nbln->shift & NBNXN_CI_HALF_LJ(0)) || !do_LJ) && do_coul;
 398
 399 #ifdef ENERGY_GROUPS
 400         egps_i = nbat->energrp[ci];
 401         {
 402             int ia, egp_ia;
 403
 404             for (ia = 0; ia < UNROLLI; ia++)
 405             {
 406                 egp_ia     = (egps_i >> (ia*egps_ishift)) & egps_imask;
 407                 vvdwtp[ia] = Vvdw + egp_ia*Vstride_i;
 408                 vctp[ia]   = Vc   + egp_ia*Vstride_i;
 409             }
 410         }
 411 #endif
 412 #if defined CALC_ENERGIES
 413 #if UNROLLJ == 4
 414         if (do_coul && l_cj[nbln->cj_ind_start].cj == ci_sh)
 415 #endif
 416 #if UNROLLJ == 2
 417         if (do_coul && l_cj[nbln->cj_ind_start].cj == (ci_sh<<1))
 418 #endif
 419 #if UNROLLJ == 8
 420         if (do_coul && l_cj[nbln->cj_ind_start].cj == (ci_sh>>1))
 421 #endif
 422         {
 423             int  ia;
 424             real Vc_sub_self;
 425
 426 #ifdef CALC_COUL_RF
 427             Vc_sub_self = 0.5*ic->c_rf;
 428 #endif
 429 #ifdef CALC_COUL_TAB
 430 #ifdef TAB_FDV0
 431             Vc_sub_self = 0.5*tab_coul_F[2];
 432 #else
 433             Vc_sub_self = 0.5*tab_coul_V[0];
 434 #endif
 435 #endif
 436 #ifdef CALC_COUL_EWALD
 437             /* beta/sqrt(pi) */
 438             Vc_sub_self = 0.5*ic->ewaldcoeff*M_2_SQRTPI;
 439 #endif
 440
 441             for (ia = 0; ia < UNROLLI; ia++)
 442             {
 443                 real qi;
 444
 445                 qi = q[sci+ia];
 446 #ifdef ENERGY_GROUPS
 447                 vctp[ia][((egps_i>>(ia*egps_ishift)) & egps_imask)*egps_jstride]
 448 #else
 449                 Vc[0]
 450 #endif
 451                     -= facel*qi*qi*Vc_sub_self;
 452             }
 453         }
 454 #endif
 455
 456         /* Load i atom data */
 457         sciy             = scix + STRIDE;
 458         sciz             = sciy + STRIDE;
 459         ix_S0            = gmx_add_pr(gmx_load1_pr(x+scix), shX_S);
 460         ix_S1            = gmx_add_pr(gmx_load1_pr(x+scix+1), shX_S);
 461         ix_S2            = gmx_add_pr(gmx_load1_pr(x+scix+2), shX_S);
 462         ix_S3            = gmx_add_pr(gmx_load1_pr(x+scix+3), shX_S);
 463         iy_S0            = gmx_add_pr(gmx_load1_pr(x+sciy), shY_S);
 464         iy_S1            = gmx_add_pr(gmx_load1_pr(x+sciy+1), shY_S);
 465         iy_S2            = gmx_add_pr(gmx_load1_pr(x+sciy+2), shY_S);
 466         iy_S3            = gmx_add_pr(gmx_load1_pr(x+sciy+3), shY_S);
 467         iz_S0            = gmx_add_pr(gmx_load1_pr(x+sciz), shZ_S);
 468         iz_S1            = gmx_add_pr(gmx_load1_pr(x+sciz+1), shZ_S);
 469         iz_S2            = gmx_add_pr(gmx_load1_pr(x+sciz+2), shZ_S);
 470         iz_S3            = gmx_add_pr(gmx_load1_pr(x+sciz+3), shZ_S);
 471
 472         if (do_coul)
 473         {
 474             iq_S0      = gmx_set1_pr(facel*q[sci]);
 475             iq_S1      = gmx_set1_pr(facel*q[sci+1]);
 476             iq_S2      = gmx_set1_pr(facel*q[sci+2]);
 477             iq_S3      = gmx_set1_pr(facel*q[sci+3]);
 478         }
 479
 480 #ifdef LJ_COMB_LB
 481         hsig_i_S0      = gmx_load1_pr(ljc+sci2+0);
 482         hsig_i_S1      = gmx_load1_pr(ljc+sci2+1);
 483         hsig_i_S2      = gmx_load1_pr(ljc+sci2+2);
 484         hsig_i_S3      = gmx_load1_pr(ljc+sci2+3);
 485         seps_i_S0      = gmx_load1_pr(ljc+sci2+STRIDE+0);
 486         seps_i_S1      = gmx_load1_pr(ljc+sci2+STRIDE+1);
 487         seps_i_S2      = gmx_load1_pr(ljc+sci2+STRIDE+2);
 488         seps_i_S3      = gmx_load1_pr(ljc+sci2+STRIDE+3);
 489 #else
 490 #ifdef LJ_COMB_GEOM
 491         c6s_S0         = gmx_load1_pr(ljc+sci2+0);
 492         c6s_S1         = gmx_load1_pr(ljc+sci2+1);
 493         if (!half_LJ)
 494         {
 495             c6s_S2     = gmx_load1_pr(ljc+sci2+2);
 496             c6s_S3     = gmx_load1_pr(ljc+sci2+3);
 497         }
 498         c12s_S0        = gmx_load1_pr(ljc+sci2+STRIDE+0);
 499         c12s_S1        = gmx_load1_pr(ljc+sci2+STRIDE+1);
 500         if (!half_LJ)
 501         {
 502             c12s_S2    = gmx_load1_pr(ljc+sci2+STRIDE+2);
 503             c12s_S3    = gmx_load1_pr(ljc+sci2+STRIDE+3);
 504         }
 505 #else
 506         nbfp0     = nbfp_ptr + type[sci  ]*nbat->ntype*nbfp_stride;
 507         nbfp1     = nbfp_ptr + type[sci+1]*nbat->ntype*nbfp_stride;
 508         if (!half_LJ)
 509         {
 510             nbfp2 = nbfp_ptr + type[sci+2]*nbat->ntype*nbfp_stride;
 511             nbfp3 = nbfp_ptr + type[sci+3]*nbat->ntype*nbfp_stride;
 512         }
 513 #endif
 514 #endif
 515
 516         /* Zero the potential energy for this list */
 517         Vvdwtot_S        = gmx_setzero_pr();
 518         vctot_S          = gmx_setzero_pr();
 519
 520         /* Clear i atom forces */
 521         fix_S0           = gmx_setzero_pr();
 522         fix_S1           = gmx_setzero_pr();
 523         fix_S2           = gmx_setzero_pr();
 524         fix_S3           = gmx_setzero_pr();
 525         fiy_S0           = gmx_setzero_pr();
 526         fiy_S1           = gmx_setzero_pr();
 527         fiy_S2           = gmx_setzero_pr();
 528         fiy_S3           = gmx_setzero_pr();
 529         fiz_S0           = gmx_setzero_pr();
 530         fiz_S1           = gmx_setzero_pr();
 531         fiz_S2           = gmx_setzero_pr();
 532         fiz_S3           = gmx_setzero_pr();
 533
 534         cjind = cjind0;
 535
 536         /* Currently all kernels use (at least half) LJ */
 537 #define CALC_LJ
 538         if (half_LJ)
 539         {
 540 #define CALC_COULOMB
 541 #define HALF_LJ
 542 #define CHECK_EXCLS
 543             while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
 544             {
 545 #include "nbnxn_kernel_simd_4xn_inner.h"
 546                 cjind++;
 547             }
 548 #undef CHECK_EXCLS
 549             for (; (cjind < cjind1); cjind++)
 550             {
 551 #include "nbnxn_kernel_simd_4xn_inner.h"
 552             }
 553 #undef HALF_LJ
 554 #undef CALC_COULOMB
 555         }
 556         else if (do_coul)
 557         {
 558 #define CALC_COULOMB
 559 #define CHECK_EXCLS
 560             while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
 561             {
 562 #include "nbnxn_kernel_simd_4xn_inner.h"
 563                 cjind++;
 564             }
 565 #undef CHECK_EXCLS
 566             for (; (cjind < cjind1); cjind++)
 567             {
 568 #include "nbnxn_kernel_simd_4xn_inner.h"
 569             }
 570 #undef CALC_COULOMB
 571         }
 572         else
 573         {
 574 #define CHECK_EXCLS
 575             while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
 576             {
 577 #include "nbnxn_kernel_simd_4xn_inner.h"
 578                 cjind++;
 579             }
 580 #undef CHECK_EXCLS
 581             for (; (cjind < cjind1); cjind++)
 582             {
 583 #include "nbnxn_kernel_simd_4xn_inner.h"
 584             }
 585         }
 586 #undef CALC_LJ
 587         ninner += cjind1 - cjind0;
 588
 589         /* Add accumulated i-forces to the force array */
 590 #if UNROLLJ >= 4
 591         fix_S = gmx_mm_transpose_sum4_pr(fix_S0, fix_S1, fix_S2, fix_S3);
 592         gmx_store_pr4(f+scix, gmx_add_pr4(fix_S, gmx_load_pr4(f+scix)));
 593
 594         fiy_S = gmx_mm_transpose_sum4_pr(fiy_S0, fiy_S1, fiy_S2, fiy_S3);
 595         gmx_store_pr4(f+sciy, gmx_add_pr4(fiy_S, gmx_load_pr4(f+sciy)));
 596
 597         fiz_S = gmx_mm_transpose_sum4_pr(fiz_S0, fiz_S1, fiz_S2, fiz_S3);
 598         gmx_store_pr4(f+sciz, gmx_add_pr4(fiz_S, gmx_load_pr4(f+sciz)));
 599
 600 #ifdef CALC_SHIFTFORCES
 601         gmx_store_pr4(shf, fix_S);
 602         fshift[ish3+0] += SUM_SIMD4(shf);
 603         gmx_store_pr4(shf, fiy_S);
 604         fshift[ish3+1] += SUM_SIMD4(shf);
 605         gmx_store_pr4(shf, fiz_S);
 606         fshift[ish3+2] += SUM_SIMD4(shf);
 607 #endif
 608 #else
 609         fix0_S = gmx_mm_transpose_sum2_pr(fix_S0, fix_S1);
 610         gmx_store_pr(f+scix, gmx_add_pr(fix0_S, gmx_load_pr(f+scix)));
 611         fix2_S = gmx_mm_transpose_sum2_pr(fix_S2, fix_S3);
 612         gmx_store_pr(f+scix+2, gmx_add_pr(fix2_S, gmx_load_pr(f+scix+2)));
 613
 614         fiy0_S = gmx_mm_transpose_sum2_pr(fiy_S0, fiy_S1);
 615         gmx_store_pr(f+sciy, gmx_add_pr(fiy0_S, gmx_load_pr(f+sciy)));
 616         fiy2_S = gmx_mm_transpose_sum2_pr(fiy_S2, fiy_S3);
 617         gmx_store_pr(f+sciy+2, gmx_add_pr(fiy2_S, gmx_load_pr(f+sciy+2)));
 618
 619         fiz0_S = gmx_mm_transpose_sum2_pr(fiz_S0, fiz_S1);
 620         gmx_store_pr(f+sciz, gmx_add_pr(fiz0_S, gmx_load_pr(f+sciz)));
 621         fiz2_S = gmx_mm_transpose_sum2_pr(fiz_S2, fiz_S3);
 622         gmx_store_pr(f+sciz+2, gmx_add_pr(fiz2_S, gmx_load_pr(f+sciz+2)));
 623
 624 #ifdef CALC_SHIFTFORCES
 625         gmx_store_pr(shf, gmx_add_pr(fix0_S, fix2_S));
 626         fshift[ish3+0] += shf[0] + shf[1];
 627         gmx_store_pr(shf, gmx_add_pr(fiy0_S, fiy2_S));
 628         fshift[ish3+1] += shf[0] + shf[1];
 629         gmx_store_pr(shf, gmx_add_pr(fiz0_S, fiz2_S));
 630         fshift[ish3+2] += shf[0] + shf[1];
 631 #endif
 632 #endif
 633
 634 #ifdef CALC_ENERGIES
 635         if (do_coul)
 636         {
 637             gmx_store_pr(tmpsum, vctot_S);
 638             *Vc += SUM_SIMD(tmpsum);
 639         }
 640
 641         gmx_store_pr(tmpsum, Vvdwtot_S);
 642         *Vvdw += SUM_SIMD(tmpsum);
 643 #endif
 644
 645         /* Outer loop uses 6 flops/iteration */
 646     }
 647
 648 #ifdef COUNT_PAIRS
 649     printf("atom pairs %d\n", npair);
 650 #endif
 651 }