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