src/gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_simd_4xn_inner.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, 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 /* This is the innermost loop contents for the 4 x N atom SIMD kernel.
  39  * This flavor of the kernel calculates interactions of 4 i-atoms
  40  * with N j-atoms stored in N wide SIMD registers.
  41  */
  42
  43
  44 /* When calculating RF or Ewald interactions we calculate the electrostatic
  45  * forces on excluded atom pairs here in the non-bonded loops.
  46  * But when energies and/or virial is required we calculate them
  47  * separately to as then it is easier to separate the energy and virial
  48  * contributions.
  49  */
  50 #if defined CHECK_EXCLS && defined CALC_COULOMB
  51 #define EXCL_FORCES
  52 #endif
  53
  54 /* Without exclusions and energies we only need to mask the cut-off,
  55  * this can be faster when we have defined gmx_blendv_pr, i.e. an instruction
  56  * that selects from two SIMD registers based on the contents of a third.
  57  */
  58 #if !(defined CHECK_EXCLS || defined CALC_ENERGIES) && defined GMX_HAVE_SIMD_BLENDV && !defined COUNT_PAIRS
  59 /* With RF and tabulated Coulomb we replace cmp+and with sub+blendv.
  60  * With gcc this is slower, except for RF on Sandy Bridge.
  61  * Tested with gcc 4.6.2, 4.6.3 and 4.7.1.
  62  */
  63 #if (defined CALC_COUL_RF || defined CALC_COUL_TAB) && (!defined __GNUC__ || (defined CALC_COUL_RF && defined GMX_X86_AVX_256))
  64 #define NBNXN_CUTOFF_USE_BLENDV
  65 #endif
  66 /* With analytical Ewald we replace cmp+and+and with sub+blendv+blendv.
  67  * This is only faster with icc on Sandy Bridge (PS kernel slower than gcc 4.7).
  68  * Tested with icc 13.
  69  */
  70 #if defined CALC_COUL_EWALD && defined __INTEL_COMPILER && defined GMX_X86_AVX_256
  71 #define NBNXN_CUTOFF_USE_BLENDV
  72 #endif
  73 #endif
  74
  75 {
  76     int        cj, aj, ajx, ajy, ajz;
  77
  78 #ifdef ENERGY_GROUPS
  79     /* Energy group indices for two atoms packed into one int */
  80     int        egp_jj[UNROLLJ/2];
  81 #endif
  82
  83 #ifdef CHECK_EXCLS
  84     /* Interaction (non-exclusion) mask of all 1's or 0's */
  85     gmx_mm_pr  int_S0;
  86     gmx_mm_pr  int_S1;
  87     gmx_mm_pr  int_S2;
  88     gmx_mm_pr  int_S3;
  89 #endif
  90
  91     gmx_mm_pr  jx_S, jy_S, jz_S;
  92     gmx_mm_pr  dx_S0, dy_S0, dz_S0;
  93     gmx_mm_pr  dx_S1, dy_S1, dz_S1;
  94     gmx_mm_pr  dx_S2, dy_S2, dz_S2;
  95     gmx_mm_pr  dx_S3, dy_S3, dz_S3;
  96     gmx_mm_pr  tx_S0, ty_S0, tz_S0;
  97     gmx_mm_pr  tx_S1, ty_S1, tz_S1;
  98     gmx_mm_pr  tx_S2, ty_S2, tz_S2;
  99     gmx_mm_pr  tx_S3, ty_S3, tz_S3;
 100     gmx_mm_pr  rsq_S0, rinv_S0, rinvsq_S0;
 101     gmx_mm_pr  rsq_S1, rinv_S1, rinvsq_S1;
 102     gmx_mm_pr  rsq_S2, rinv_S2, rinvsq_S2;
 103     gmx_mm_pr  rsq_S3, rinv_S3, rinvsq_S3;
 104 #ifndef NBNXN_CUTOFF_USE_BLENDV
 105     /* wco: within cut-off, mask of all 1's or 0's */
 106     gmx_mm_pr  wco_S0;
 107     gmx_mm_pr  wco_S1;
 108     gmx_mm_pr  wco_S2;
 109     gmx_mm_pr  wco_S3;
 110 #endif
 111 #ifdef VDW_CUTOFF_CHECK
 112     gmx_mm_pr  wco_vdw_S0;
 113     gmx_mm_pr  wco_vdw_S1;
 114 #ifndef HALF_LJ
 115     gmx_mm_pr  wco_vdw_S2;
 116     gmx_mm_pr  wco_vdw_S3;
 117 #endif
 118 #endif
 119 #ifdef CALC_COULOMB
 120 #ifdef CHECK_EXCLS
 121     /* 1/r masked with the interaction mask */
 122     gmx_mm_pr  rinv_ex_S0;
 123     gmx_mm_pr  rinv_ex_S1;
 124     gmx_mm_pr  rinv_ex_S2;
 125     gmx_mm_pr  rinv_ex_S3;
 126 #endif
 127     gmx_mm_pr  jq_S;
 128     gmx_mm_pr  qq_S0;
 129     gmx_mm_pr  qq_S1;
 130     gmx_mm_pr  qq_S2;
 131     gmx_mm_pr  qq_S3;
 132 #ifdef CALC_COUL_TAB
 133     /* The force (PME mesh force) we need to subtract from 1/r^2 */
 134     gmx_mm_pr  fsub_S0;
 135     gmx_mm_pr  fsub_S1;
 136     gmx_mm_pr  fsub_S2;
 137     gmx_mm_pr  fsub_S3;
 138 #endif
 139 #ifdef CALC_COUL_EWALD
 140     gmx_mm_pr  brsq_S0, brsq_S1, brsq_S2, brsq_S3;
 141     gmx_mm_pr  ewcorr_S0, ewcorr_S1, ewcorr_S2, ewcorr_S3;
 142 #endif
 143
 144     /* frcoul = (1/r - fsub)*r */
 145     gmx_mm_pr  frcoul_S0;
 146     gmx_mm_pr  frcoul_S1;
 147     gmx_mm_pr  frcoul_S2;
 148     gmx_mm_pr  frcoul_S3;
 149 #ifdef CALC_COUL_TAB
 150     /* For tables: r, rs=r/sp, rf=floor(rs), frac=rs-rf */
 151     gmx_mm_pr  r_S0, rs_S0, rf_S0, frac_S0;
 152     gmx_mm_pr  r_S1, rs_S1, rf_S1, frac_S1;
 153     gmx_mm_pr  r_S2, rs_S2, rf_S2, frac_S2;
 154     gmx_mm_pr  r_S3, rs_S3, rf_S3, frac_S3;
 155     /* Table index: rs truncated to an int */
 156     gmx_epi32  ti_S0, ti_S1, ti_S2, ti_S3;
 157     /* Linear force table values */
 158     gmx_mm_pr  ctab0_S0, ctab1_S0;
 159     gmx_mm_pr  ctab0_S1, ctab1_S1;
 160     gmx_mm_pr  ctab0_S2, ctab1_S2;
 161     gmx_mm_pr  ctab0_S3, ctab1_S3;
 162 #ifdef CALC_ENERGIES
 163     /* Quadratic energy table value */
 164     gmx_mm_pr  ctabv_S0;
 165     gmx_mm_pr  ctabv_S1;
 166     gmx_mm_pr  ctabv_S2;
 167     gmx_mm_pr  ctabv_S3;
 168 #endif
 169 #endif
 170 #if defined CALC_ENERGIES && (defined CALC_COUL_EWALD || defined CALC_COUL_TAB)
 171     /* The potential (PME mesh) we need to subtract from 1/r */
 172     gmx_mm_pr  vc_sub_S0;
 173     gmx_mm_pr  vc_sub_S1;
 174     gmx_mm_pr  vc_sub_S2;
 175     gmx_mm_pr  vc_sub_S3;
 176 #endif
 177 #ifdef CALC_ENERGIES
 178     /* Electrostatic potential */
 179     gmx_mm_pr  vcoul_S0;
 180     gmx_mm_pr  vcoul_S1;
 181     gmx_mm_pr  vcoul_S2;
 182     gmx_mm_pr  vcoul_S3;
 183 #endif
 184 #endif
 185     /* The force times 1/r */
 186     gmx_mm_pr  fscal_S0;
 187     gmx_mm_pr  fscal_S1;
 188     gmx_mm_pr  fscal_S2;
 189     gmx_mm_pr  fscal_S3;
 190
 191 #ifdef CALC_LJ
 192 #ifdef LJ_COMB_LB
 193     /* LJ sigma_j/2 and sqrt(epsilon_j) */
 194     gmx_mm_pr  hsig_j_S, seps_j_S;
 195     /* LJ sigma_ij and epsilon_ij */
 196     gmx_mm_pr  sig_S0, eps_S0;
 197     gmx_mm_pr  sig_S1, eps_S1;
 198 #ifndef HALF_LJ
 199     gmx_mm_pr  sig_S2, eps_S2;
 200     gmx_mm_pr  sig_S3, eps_S3;
 201 #endif
 202 #ifdef CALC_ENERGIES
 203     gmx_mm_pr  sig2_S0, sig6_S0;
 204     gmx_mm_pr  sig2_S1, sig6_S1;
 205 #ifndef HALF_LJ
 206     gmx_mm_pr  sig2_S2, sig6_S2;
 207     gmx_mm_pr  sig2_S3, sig6_S3;
 208 #endif
 209 #endif /* LJ_COMB_LB */
 210 #endif /* CALC_LJ */
 211
 212 #ifdef LJ_COMB_GEOM
 213     gmx_mm_pr  c6s_j_S, c12s_j_S;
 214 #endif
 215
 216 #if defined LJ_COMB_GEOM || defined LJ_COMB_LB
 217     /* Index for loading LJ parameters, complicated when interleaving */
 218     int         aj2;
 219 #endif
 220
 221 #ifndef FIX_LJ_C
 222     /* LJ C6 and C12 parameters, used with geometric comb. rule */
 223     gmx_mm_pr  c6_S0, c12_S0;
 224     gmx_mm_pr  c6_S1, c12_S1;
 225 #ifndef HALF_LJ
 226     gmx_mm_pr  c6_S2, c12_S2;
 227     gmx_mm_pr  c6_S3, c12_S3;
 228 #endif
 229 #endif
 230
 231     /* Intermediate variables for LJ calculation */
 232 #ifndef LJ_COMB_LB
 233     gmx_mm_pr  rinvsix_S0;
 234     gmx_mm_pr  rinvsix_S1;
 235 #ifndef HALF_LJ
 236     gmx_mm_pr  rinvsix_S2;
 237     gmx_mm_pr  rinvsix_S3;
 238 #endif
 239 #endif
 240 #ifdef LJ_COMB_LB
 241     gmx_mm_pr  sir_S0, sir2_S0, sir6_S0;
 242     gmx_mm_pr  sir_S1, sir2_S1, sir6_S1;
 243 #ifndef HALF_LJ
 244     gmx_mm_pr  sir_S2, sir2_S2, sir6_S2;
 245     gmx_mm_pr  sir_S3, sir2_S3, sir6_S3;
 246 #endif
 247 #endif
 248
 249     gmx_mm_pr  FrLJ6_S0, FrLJ12_S0;
 250     gmx_mm_pr  FrLJ6_S1, FrLJ12_S1;
 251 #ifndef HALF_LJ
 252     gmx_mm_pr  FrLJ6_S2, FrLJ12_S2;
 253     gmx_mm_pr  FrLJ6_S3, FrLJ12_S3;
 254 #endif
 255 #ifdef CALC_ENERGIES
 256     gmx_mm_pr  VLJ6_S0, VLJ12_S0, VLJ_S0;
 257     gmx_mm_pr  VLJ6_S1, VLJ12_S1, VLJ_S1;
 258 #ifndef HALF_LJ
 259     gmx_mm_pr  VLJ6_S2, VLJ12_S2, VLJ_S2;
 260     gmx_mm_pr  VLJ6_S3, VLJ12_S3, VLJ_S3;
 261 #endif
 262 #endif
 263 #endif /* CALC_LJ */
 264
 265     /* j-cluster index */
 266     cj            = l_cj[cjind].cj;
 267
 268     /* Atom indices (of the first atom in the cluster) */
 269     aj            = cj*UNROLLJ;
 270 #if defined CALC_LJ && (defined LJ_COMB_GEOM || defined LJ_COMB_LB)
 271 #if UNROLLJ == STRIDE
 272     aj2           = aj*2;
 273 #else
 274     aj2           = (cj>>1)*2*STRIDE + (cj & 1)*UNROLLJ;
 275 #endif
 276 #endif
 277 #if UNROLLJ == STRIDE
 278     ajx           = aj*DIM;
 279 #else
 280     ajx           = (cj>>1)*DIM*STRIDE + (cj & 1)*UNROLLJ;
 281 #endif
 282     ajy           = ajx + STRIDE;
 283     ajz           = ajy + STRIDE;
 284
 285 #ifdef CHECK_EXCLS
 286 #ifdef gmx_checkbitmask_epi32
 287     {
 288         /* Integer mask set and operations, cast result to real */
 289         gmx_epi32 mask_pr_S = gmx_set1_epi32(l_cj[cjind].excl);
 290
 291         int_S0  = gmx_castsi_pr(gmx_checkbitmask_epi32(mask_pr_S, mask_S0));
 292         int_S1  = gmx_castsi_pr(gmx_checkbitmask_epi32(mask_pr_S, mask_S1));
 293         int_S2  = gmx_castsi_pr(gmx_checkbitmask_epi32(mask_pr_S, mask_S2));
 294         int_S3  = gmx_castsi_pr(gmx_checkbitmask_epi32(mask_pr_S, mask_S3));
 295     }
 296 #else
 297     {
 298         /* Integer mask set, cast to real and real mask operations */
 299         gmx_mm_pr mask_pr_S = gmx_castsi_pr(gmx_set1_epi32(l_cj[cjind].excl));
 300
 301         int_S0  = gmx_checkbitmask_pr(mask_pr_S, mask_S0);
 302         int_S1  = gmx_checkbitmask_pr(mask_pr_S, mask_S1);
 303         int_S2  = gmx_checkbitmask_pr(mask_pr_S, mask_S2);
 304         int_S3  = gmx_checkbitmask_pr(mask_pr_S, mask_S3);
 305     }
 306 #endif
 307 #endif
 308
 309     /* load j atom coordinates */
 310     jx_S        = gmx_load_pr(x+ajx);
 311     jy_S        = gmx_load_pr(x+ajy);
 312     jz_S        = gmx_load_pr(x+ajz);
 313
 314     /* Calculate distance */
 315     dx_S0       = gmx_sub_pr(ix_S0, jx_S);
 316     dy_S0       = gmx_sub_pr(iy_S0, jy_S);
 317     dz_S0       = gmx_sub_pr(iz_S0, jz_S);
 318     dx_S1       = gmx_sub_pr(ix_S1, jx_S);
 319     dy_S1       = gmx_sub_pr(iy_S1, jy_S);
 320     dz_S1       = gmx_sub_pr(iz_S1, jz_S);
 321     dx_S2       = gmx_sub_pr(ix_S2, jx_S);
 322     dy_S2       = gmx_sub_pr(iy_S2, jy_S);
 323     dz_S2       = gmx_sub_pr(iz_S2, jz_S);
 324     dx_S3       = gmx_sub_pr(ix_S3, jx_S);
 325     dy_S3       = gmx_sub_pr(iy_S3, jy_S);
 326     dz_S3       = gmx_sub_pr(iz_S3, jz_S);
 327
 328     /* rsq = dx*dx+dy*dy+dz*dz */
 329     rsq_S0      = gmx_calc_rsq_pr(dx_S0, dy_S0, dz_S0);
 330     rsq_S1      = gmx_calc_rsq_pr(dx_S1, dy_S1, dz_S1);
 331     rsq_S2      = gmx_calc_rsq_pr(dx_S2, dy_S2, dz_S2);
 332     rsq_S3      = gmx_calc_rsq_pr(dx_S3, dy_S3, dz_S3);
 333
 334 #ifndef NBNXN_CUTOFF_USE_BLENDV
 335     wco_S0      = gmx_cmplt_pr(rsq_S0, rc2_S);
 336     wco_S1      = gmx_cmplt_pr(rsq_S1, rc2_S);
 337     wco_S2      = gmx_cmplt_pr(rsq_S2, rc2_S);
 338     wco_S3      = gmx_cmplt_pr(rsq_S3, rc2_S);
 339 #endif
 340
 341 #ifdef CHECK_EXCLS
 342 #ifdef EXCL_FORCES
 343     /* Only remove the (sub-)diagonal to avoid double counting */
 344 #if UNROLLJ == UNROLLI
 345     if (cj == ci_sh)
 346     {
 347         wco_S0  = gmx_and_pr(wco_S0, diag_S0);
 348         wco_S1  = gmx_and_pr(wco_S1, diag_S1);
 349         wco_S2  = gmx_and_pr(wco_S2, diag_S2);
 350         wco_S3  = gmx_and_pr(wco_S3, diag_S3);
 351     }
 352 #else
 353 #if UNROLLJ < UNROLLI
 354     if (cj == ci_sh*2)
 355     {
 356         wco_S0  = gmx_and_pr(wco_S0, diag0_S0);
 357         wco_S1  = gmx_and_pr(wco_S1, diag0_S1);
 358         wco_S2  = gmx_and_pr(wco_S2, diag0_S2);
 359         wco_S3  = gmx_and_pr(wco_S3, diag0_S3);
 360     }
 361     if (cj == ci_sh*2 + 1)
 362     {
 363         wco_S0  = gmx_and_pr(wco_S0, diag1_S0);
 364         wco_S1  = gmx_and_pr(wco_S1, diag1_S1);
 365         wco_S2  = gmx_and_pr(wco_S2, diag1_S2);
 366         wco_S3  = gmx_and_pr(wco_S3, diag1_S3);
 367     }
 368 #else
 369     if (cj*2 == ci_sh)
 370     {
 371         wco_S0  = gmx_and_pr(wco_S0, diag0_S0);
 372         wco_S1  = gmx_and_pr(wco_S1, diag0_S1);
 373         wco_S2  = gmx_and_pr(wco_S2, diag0_S2);
 374         wco_S3  = gmx_and_pr(wco_S3, diag0_S3);
 375     }
 376     else if (cj*2 + 1 == ci_sh)
 377     {
 378         wco_S0  = gmx_and_pr(wco_S0, diag1_S0);
 379         wco_S1  = gmx_and_pr(wco_S1, diag1_S1);
 380         wco_S2  = gmx_and_pr(wco_S2, diag1_S2);
 381         wco_S3  = gmx_and_pr(wco_S3, diag1_S3);
 382     }
 383 #endif
 384 #endif
 385 #else /* EXCL_FORCES */
 386     /* No exclusion forces: remove all excluded atom pairs from the list */
 387     wco_S0      = gmx_and_pr(wco_S0, int_S0);
 388     wco_S1      = gmx_and_pr(wco_S1, int_S1);
 389     wco_S2      = gmx_and_pr(wco_S2, int_S2);
 390     wco_S3      = gmx_and_pr(wco_S3, int_S3);
 391 #endif
 392 #endif
 393
 394 #ifdef COUNT_PAIRS
 395     {
 396         int  i, j;
 397         real tmpa[2*GMX_SIMD_WIDTH_HERE], *tmp;
 398         tmp = gmx_simd_align_real(tmpa);
 399         for (i = 0; i < UNROLLI; i++)
 400         {
 401             gmx_store_pr(tmp, i == 0 ? wco_S0 : (i == 1 ? wco_S1 : (i == 2 ? wco_S2 : wco_S3)));
 402             for (j = 0; j < UNROLLJ; j++)
 403             {
 404                 if (!(tmp[j] == 0))
 405                 {
 406                     npair++;
 407                 }
 408             }
 409         }
 410     }
 411 #endif
 412
 413 #ifdef CHECK_EXCLS
 414     /* For excluded pairs add a small number to avoid r^-6 = NaN */
 415     rsq_S0      = gmx_add_pr(rsq_S0, gmx_andnot_pr(int_S0, avoid_sing_S));
 416     rsq_S1      = gmx_add_pr(rsq_S1, gmx_andnot_pr(int_S1, avoid_sing_S));
 417     rsq_S2      = gmx_add_pr(rsq_S2, gmx_andnot_pr(int_S2, avoid_sing_S));
 418     rsq_S3      = gmx_add_pr(rsq_S3, gmx_andnot_pr(int_S3, avoid_sing_S));
 419 #endif
 420
 421     /* Calculate 1/r */
 422 #ifndef GMX_DOUBLE
 423     rinv_S0     = gmx_invsqrt_pr(rsq_S0);
 424     rinv_S1     = gmx_invsqrt_pr(rsq_S1);
 425     rinv_S2     = gmx_invsqrt_pr(rsq_S2);
 426     rinv_S3     = gmx_invsqrt_pr(rsq_S3);
 427 #else
 428     GMX_MM_INVSQRT2_PD(rsq_S0, rsq_S1, rinv_S0, rinv_S1);
 429     GMX_MM_INVSQRT2_PD(rsq_S2, rsq_S3, rinv_S2, rinv_S3);
 430 #endif
 431
 432 #ifdef CALC_COULOMB
 433     /* Load parameters for j atom */
 434     jq_S        = gmx_load_pr(q+aj);
 435     qq_S0       = gmx_mul_pr(iq_S0, jq_S);
 436     qq_S1       = gmx_mul_pr(iq_S1, jq_S);
 437     qq_S2       = gmx_mul_pr(iq_S2, jq_S);
 438     qq_S3       = gmx_mul_pr(iq_S3, jq_S);
 439 #endif
 440
 441 #ifdef CALC_LJ
 442
 443 #if !defined LJ_COMB_GEOM && !defined LJ_COMB_LB && !defined FIX_LJ_C
 444     load_lj_pair_params(nbfp0, type, aj, c6_S0, c12_S0);
 445     load_lj_pair_params(nbfp1, type, aj, c6_S1, c12_S1);
 446 #ifndef HALF_LJ
 447     load_lj_pair_params(nbfp2, type, aj, c6_S2, c12_S2);
 448     load_lj_pair_params(nbfp3, type, aj, c6_S3, c12_S3);
 449 #endif
 450 #endif /* not defined any LJ rule */
 451
 452 #ifdef LJ_COMB_GEOM
 453     c6s_j_S     = gmx_load_pr(ljc+aj2+0);
 454     c12s_j_S    = gmx_load_pr(ljc+aj2+STRIDE);
 455     c6_S0       = gmx_mul_pr(c6s_S0, c6s_j_S );
 456     c6_S1       = gmx_mul_pr(c6s_S1, c6s_j_S );
 457 #ifndef HALF_LJ
 458     c6_S2       = gmx_mul_pr(c6s_S2, c6s_j_S );
 459     c6_S3       = gmx_mul_pr(c6s_S3, c6s_j_S );
 460 #endif
 461     c12_S0      = gmx_mul_pr(c12s_S0, c12s_j_S);
 462     c12_S1      = gmx_mul_pr(c12s_S1, c12s_j_S);
 463 #ifndef HALF_LJ
 464     c12_S2      = gmx_mul_pr(c12s_S2, c12s_j_S);
 465     c12_S3      = gmx_mul_pr(c12s_S3, c12s_j_S);
 466 #endif
 467 #endif /* LJ_COMB_GEOM */
 468
 469 #ifdef LJ_COMB_LB
 470     hsig_j_S    = gmx_load_pr(ljc+aj2+0);
 471     seps_j_S    = gmx_load_pr(ljc+aj2+STRIDE);
 472
 473     sig_S0      = gmx_add_pr(hsig_i_S0, hsig_j_S);
 474     sig_S1      = gmx_add_pr(hsig_i_S1, hsig_j_S);
 475     eps_S0      = gmx_mul_pr(seps_i_S0, seps_j_S);
 476     eps_S1      = gmx_mul_pr(seps_i_S1, seps_j_S);
 477 #ifndef HALF_LJ
 478     sig_S2      = gmx_add_pr(hsig_i_S2, hsig_j_S);
 479     sig_S3      = gmx_add_pr(hsig_i_S3, hsig_j_S);
 480     eps_S2      = gmx_mul_pr(seps_i_S2, seps_j_S);
 481     eps_S3      = gmx_mul_pr(seps_i_S3, seps_j_S);
 482 #endif
 483 #endif /* LJ_COMB_LB */
 484
 485 #endif /* CALC_LJ */
 486
 487 #ifndef NBNXN_CUTOFF_USE_BLENDV
 488     rinv_S0     = gmx_blendzero_pr(rinv_S0, wco_S0);
 489     rinv_S1     = gmx_blendzero_pr(rinv_S1, wco_S1);
 490     rinv_S2     = gmx_blendzero_pr(rinv_S2, wco_S2);
 491     rinv_S3     = gmx_blendzero_pr(rinv_S3, wco_S3);
 492 #else
 493     /* We only need to mask for the cut-off: blendv is faster */
 494     rinv_S0     = gmx_blendv_pr(rinv_S0, zero_S, gmx_sub_pr(rc2_S, rsq_S0));
 495     rinv_S1     = gmx_blendv_pr(rinv_S1, zero_S, gmx_sub_pr(rc2_S, rsq_S1));
 496     rinv_S2     = gmx_blendv_pr(rinv_S2, zero_S, gmx_sub_pr(rc2_S, rsq_S2));
 497     rinv_S3     = gmx_blendv_pr(rinv_S3, zero_S, gmx_sub_pr(rc2_S, rsq_S3));
 498 #endif
 499
 500     rinvsq_S0   = gmx_mul_pr(rinv_S0, rinv_S0);
 501     rinvsq_S1   = gmx_mul_pr(rinv_S1, rinv_S1);
 502     rinvsq_S2   = gmx_mul_pr(rinv_S2, rinv_S2);
 503     rinvsq_S3   = gmx_mul_pr(rinv_S3, rinv_S3);
 504
 505 #ifdef CALC_COULOMB
 506     /* Note that here we calculate force*r, not the usual force/r.
 507      * This allows avoiding masking the reaction-field contribution,
 508      * as frcoul is later multiplied by rinvsq which has been
 509      * masked with the cut-off check.
 510      */
 511
 512 #ifdef EXCL_FORCES
 513     /* Only add 1/r for non-excluded atom pairs */
 514     rinv_ex_S0  = gmx_blendzero_pr(rinv_S0, int_S0);
 515     rinv_ex_S1  = gmx_blendzero_pr(rinv_S1, int_S1);
 516     rinv_ex_S2  = gmx_blendzero_pr(rinv_S2, int_S2);
 517     rinv_ex_S3  = gmx_blendzero_pr(rinv_S3, int_S3);
 518 #else
 519     /* No exclusion forces, we always need 1/r */
 520 #define     rinv_ex_S0    rinv_S0
 521 #define     rinv_ex_S1    rinv_S1
 522 #define     rinv_ex_S2    rinv_S2
 523 #define     rinv_ex_S3    rinv_S3
 524 #endif
 525
 526 #ifdef CALC_COUL_RF
 527     /* Electrostatic interactions */
 528     frcoul_S0   = gmx_mul_pr(qq_S0, gmx_add_pr(rinv_ex_S0, gmx_mul_pr(rsq_S0, mrc_3_S)));
 529     frcoul_S1   = gmx_mul_pr(qq_S1, gmx_add_pr(rinv_ex_S1, gmx_mul_pr(rsq_S1, mrc_3_S)));
 530     frcoul_S2   = gmx_mul_pr(qq_S2, gmx_add_pr(rinv_ex_S2, gmx_mul_pr(rsq_S2, mrc_3_S)));
 531     frcoul_S3   = gmx_mul_pr(qq_S3, gmx_add_pr(rinv_ex_S3, gmx_mul_pr(rsq_S3, mrc_3_S)));
 532
 533 #ifdef CALC_ENERGIES
 534     vcoul_S0    = gmx_mul_pr(qq_S0, gmx_add_pr(rinv_ex_S0, gmx_add_pr(gmx_mul_pr(rsq_S0, hrc_3_S), moh_rc_S)));
 535     vcoul_S1    = gmx_mul_pr(qq_S1, gmx_add_pr(rinv_ex_S1, gmx_add_pr(gmx_mul_pr(rsq_S1, hrc_3_S), moh_rc_S)));
 536     vcoul_S2    = gmx_mul_pr(qq_S2, gmx_add_pr(rinv_ex_S2, gmx_add_pr(gmx_mul_pr(rsq_S2, hrc_3_S), moh_rc_S)));
 537     vcoul_S3    = gmx_mul_pr(qq_S3, gmx_add_pr(rinv_ex_S3, gmx_add_pr(gmx_mul_pr(rsq_S3, hrc_3_S), moh_rc_S)));
 538 #endif
 539 #endif
 540
 541 #ifdef CALC_COUL_EWALD
 542     /* We need to mask (or limit) rsq for the cut-off,
 543      * as large distances can cause an overflow in gmx_pmecorrF/V.
 544      */
 545 #ifndef NBNXN_CUTOFF_USE_BLENDV
 546     brsq_S0     = gmx_mul_pr(beta2_S, gmx_blendzero_pr(rsq_S0, wco_S0));
 547     brsq_S1     = gmx_mul_pr(beta2_S, gmx_blendzero_pr(rsq_S1, wco_S1));
 548     brsq_S2     = gmx_mul_pr(beta2_S, gmx_blendzero_pr(rsq_S2, wco_S2));
 549     brsq_S3     = gmx_mul_pr(beta2_S, gmx_blendzero_pr(rsq_S3, wco_S3));
 550 #else
 551     /* Strangely, putting mul on a separate line is slower (icc 13) */
 552     brsq_S0     = gmx_mul_pr(beta2_S, gmx_blendv_pr(rsq_S0, zero_S, gmx_sub_pr(rc2_S, rsq_S0)));
 553     brsq_S1     = gmx_mul_pr(beta2_S, gmx_blendv_pr(rsq_S1, zero_S, gmx_sub_pr(rc2_S, rsq_S1)));
 554     brsq_S2     = gmx_mul_pr(beta2_S, gmx_blendv_pr(rsq_S2, zero_S, gmx_sub_pr(rc2_S, rsq_S2)));
 555     brsq_S3     = gmx_mul_pr(beta2_S, gmx_blendv_pr(rsq_S3, zero_S, gmx_sub_pr(rc2_S, rsq_S3)));
 556 #endif
 557     ewcorr_S0   = gmx_mul_pr(gmx_pmecorrF_pr(brsq_S0), beta_S);
 558     ewcorr_S1   = gmx_mul_pr(gmx_pmecorrF_pr(brsq_S1), beta_S);
 559     ewcorr_S2   = gmx_mul_pr(gmx_pmecorrF_pr(brsq_S2), beta_S);
 560     ewcorr_S3   = gmx_mul_pr(gmx_pmecorrF_pr(brsq_S3), beta_S);
 561     frcoul_S0   = gmx_mul_pr(qq_S0, gmx_add_pr(rinv_ex_S0, gmx_mul_pr(ewcorr_S0, brsq_S0)));
 562     frcoul_S1   = gmx_mul_pr(qq_S1, gmx_add_pr(rinv_ex_S1, gmx_mul_pr(ewcorr_S1, brsq_S1)));
 563     frcoul_S2   = gmx_mul_pr(qq_S2, gmx_add_pr(rinv_ex_S2, gmx_mul_pr(ewcorr_S2, brsq_S2)));
 564     frcoul_S3   = gmx_mul_pr(qq_S3, gmx_add_pr(rinv_ex_S3, gmx_mul_pr(ewcorr_S3, brsq_S3)));
 565
 566 #ifdef CALC_ENERGIES
 567     vc_sub_S0   = gmx_mul_pr(gmx_pmecorrV_pr(brsq_S0), beta_S);
 568     vc_sub_S1   = gmx_mul_pr(gmx_pmecorrV_pr(brsq_S1), beta_S);
 569     vc_sub_S2   = gmx_mul_pr(gmx_pmecorrV_pr(brsq_S2), beta_S);
 570     vc_sub_S3   = gmx_mul_pr(gmx_pmecorrV_pr(brsq_S3), beta_S);
 571 #endif
 572
 573 #endif /* CALC_COUL_EWALD */
 574
 575 #ifdef CALC_COUL_TAB
 576     /* Electrostatic interactions */
 577     r_S0        = gmx_mul_pr(rsq_S0, rinv_S0);
 578     r_S1        = gmx_mul_pr(rsq_S1, rinv_S1);
 579     r_S2        = gmx_mul_pr(rsq_S2, rinv_S2);
 580     r_S3        = gmx_mul_pr(rsq_S3, rinv_S3);
 581     /* Convert r to scaled table units */
 582     rs_S0       = gmx_mul_pr(r_S0, invtsp_S);
 583     rs_S1       = gmx_mul_pr(r_S1, invtsp_S);
 584     rs_S2       = gmx_mul_pr(r_S2, invtsp_S);
 585     rs_S3       = gmx_mul_pr(r_S3, invtsp_S);
 586     /* Truncate scaled r to an int */
 587     ti_S0       = gmx_cvttpr_epi32(rs_S0);
 588     ti_S1       = gmx_cvttpr_epi32(rs_S1);
 589     ti_S2       = gmx_cvttpr_epi32(rs_S2);
 590     ti_S3       = gmx_cvttpr_epi32(rs_S3);
 591 #ifdef GMX_HAVE_SIMD_FLOOR
 592     /* SSE4.1 floor is faster than gmx_cvtepi32_ps int->float cast */
 593     rf_S0       = gmx_floor_pr(rs_S0);
 594     rf_S1       = gmx_floor_pr(rs_S1);
 595     rf_S2       = gmx_floor_pr(rs_S2);
 596     rf_S3       = gmx_floor_pr(rs_S3);
 597 #else
 598     rf_S0       = gmx_cvtepi32_pr(ti_S0);
 599     rf_S1       = gmx_cvtepi32_pr(ti_S1);
 600     rf_S2       = gmx_cvtepi32_pr(ti_S2);
 601     rf_S3       = gmx_cvtepi32_pr(ti_S3);
 602 #endif
 603     frac_S0     = gmx_sub_pr(rs_S0, rf_S0);
 604     frac_S1     = gmx_sub_pr(rs_S1, rf_S1);
 605     frac_S2     = gmx_sub_pr(rs_S2, rf_S2);
 606     frac_S3     = gmx_sub_pr(rs_S3, rf_S3);
 607
 608     /* Load and interpolate table forces and possibly energies.
 609      * Force and energy can be combined in one table, stride 4: FDV0
 610      * or in two separate tables with stride 1: F and V
 611      * Currently single precision uses FDV0, double F and V.
 612      */
 613 #ifndef CALC_ENERGIES
 614     load_table_f(tab_coul_F, ti_S0, ti0, ctab0_S0, ctab1_S0);
 615     load_table_f(tab_coul_F, ti_S1, ti1, ctab0_S1, ctab1_S1);
 616     load_table_f(tab_coul_F, ti_S2, ti2, ctab0_S2, ctab1_S2);
 617     load_table_f(tab_coul_F, ti_S3, ti3, ctab0_S3, ctab1_S3);
 618 #else
 619 #ifdef TAB_FDV0
 620     load_table_f_v(tab_coul_F, ti_S0, ti0, ctab0_S0, ctab1_S0, ctabv_S0);
 621     load_table_f_v(tab_coul_F, ti_S1, ti1, ctab0_S1, ctab1_S1, ctabv_S1);
 622     load_table_f_v(tab_coul_F, ti_S2, ti2, ctab0_S2, ctab1_S2, ctabv_S2);
 623     load_table_f_v(tab_coul_F, ti_S3, ti3, ctab0_S3, ctab1_S3, ctabv_S3);
 624 #else
 625     load_table_f_v(tab_coul_F, tab_coul_V, ti_S0, ti0, ctab0_S0, ctab1_S0, ctabv_S0);
 626     load_table_f_v(tab_coul_F, tab_coul_V, ti_S1, ti1, ctab0_S1, ctab1_S1, ctabv_S1);
 627     load_table_f_v(tab_coul_F, tab_coul_V, ti_S2, ti2, ctab0_S2, ctab1_S2, ctabv_S2);
 628     load_table_f_v(tab_coul_F, tab_coul_V, ti_S3, ti3, ctab0_S3, ctab1_S3, ctabv_S3);
 629 #endif
 630 #endif
 631     fsub_S0     = gmx_add_pr(ctab0_S0, gmx_mul_pr(frac_S0, ctab1_S0));
 632     fsub_S1     = gmx_add_pr(ctab0_S1, gmx_mul_pr(frac_S1, ctab1_S1));
 633     fsub_S2     = gmx_add_pr(ctab0_S2, gmx_mul_pr(frac_S2, ctab1_S2));
 634     fsub_S3     = gmx_add_pr(ctab0_S3, gmx_mul_pr(frac_S3, ctab1_S3));
 635     frcoul_S0   = gmx_mul_pr(qq_S0, gmx_sub_pr(rinv_ex_S0, gmx_mul_pr(fsub_S0, r_S0)));
 636     frcoul_S1   = gmx_mul_pr(qq_S1, gmx_sub_pr(rinv_ex_S1, gmx_mul_pr(fsub_S1, r_S1)));
 637     frcoul_S2   = gmx_mul_pr(qq_S2, gmx_sub_pr(rinv_ex_S2, gmx_mul_pr(fsub_S2, r_S2)));
 638     frcoul_S3   = gmx_mul_pr(qq_S3, gmx_sub_pr(rinv_ex_S3, gmx_mul_pr(fsub_S3, r_S3)));
 639
 640 #ifdef CALC_ENERGIES
 641     vc_sub_S0   = gmx_add_pr(ctabv_S0, gmx_mul_pr(gmx_mul_pr(mhalfsp_S, frac_S0), gmx_add_pr(ctab0_S0, fsub_S0)));
 642     vc_sub_S1   = gmx_add_pr(ctabv_S1, gmx_mul_pr(gmx_mul_pr(mhalfsp_S, frac_S1), gmx_add_pr(ctab0_S1, fsub_S1)));
 643     vc_sub_S2   = gmx_add_pr(ctabv_S2, gmx_mul_pr(gmx_mul_pr(mhalfsp_S, frac_S2), gmx_add_pr(ctab0_S2, fsub_S2)));
 644     vc_sub_S3   = gmx_add_pr(ctabv_S3, gmx_mul_pr(gmx_mul_pr(mhalfsp_S, frac_S3), gmx_add_pr(ctab0_S3, fsub_S3)));
 645 #endif
 646 #endif /* CALC_COUL_TAB */
 647
 648 #if defined CALC_ENERGIES && (defined CALC_COUL_EWALD || defined CALC_COUL_TAB)
 649 #ifndef NO_SHIFT_EWALD
 650     /* Add Ewald potential shift to vc_sub for convenience */
 651 #ifdef CHECK_EXCLS
 652     vc_sub_S0   = gmx_add_pr(vc_sub_S0, gmx_blendzero_pr(sh_ewald_S, int_S0));
 653     vc_sub_S1   = gmx_add_pr(vc_sub_S1, gmx_blendzero_pr(sh_ewald_S, int_S1));
 654     vc_sub_S2   = gmx_add_pr(vc_sub_S2, gmx_blendzero_pr(sh_ewald_S, int_S2));
 655     vc_sub_S3   = gmx_add_pr(vc_sub_S3, gmx_blendzero_pr(sh_ewald_S, int_S3));
 656 #else
 657     vc_sub_S0   = gmx_add_pr(vc_sub_S0, sh_ewald_S);
 658     vc_sub_S1   = gmx_add_pr(vc_sub_S1, sh_ewald_S);
 659     vc_sub_S2   = gmx_add_pr(vc_sub_S2, sh_ewald_S);
 660     vc_sub_S3   = gmx_add_pr(vc_sub_S3, sh_ewald_S);
 661 #endif
 662 #endif
 663
 664     vcoul_S0    = gmx_mul_pr(qq_S0, gmx_sub_pr(rinv_ex_S0, vc_sub_S0));
 665     vcoul_S1    = gmx_mul_pr(qq_S1, gmx_sub_pr(rinv_ex_S1, vc_sub_S1));
 666     vcoul_S2    = gmx_mul_pr(qq_S2, gmx_sub_pr(rinv_ex_S2, vc_sub_S2));
 667     vcoul_S3    = gmx_mul_pr(qq_S3, gmx_sub_pr(rinv_ex_S3, vc_sub_S3));
 668
 669 #endif
 670
 671 #ifdef CALC_ENERGIES
 672     /* Mask energy for cut-off and diagonal */
 673     vcoul_S0    = gmx_blendzero_pr(vcoul_S0, wco_S0);
 674     vcoul_S1    = gmx_blendzero_pr(vcoul_S1, wco_S1);
 675     vcoul_S2    = gmx_blendzero_pr(vcoul_S2, wco_S2);
 676     vcoul_S3    = gmx_blendzero_pr(vcoul_S3, wco_S3);
 677 #endif
 678
 679 #endif /* CALC_COULOMB */
 680
 681 #ifdef CALC_LJ
 682     /* Lennard-Jones interaction */
 683
 684 #ifdef VDW_CUTOFF_CHECK
 685     wco_vdw_S0  = gmx_cmplt_pr(rsq_S0, rcvdw2_S);
 686     wco_vdw_S1  = gmx_cmplt_pr(rsq_S1, rcvdw2_S);
 687 #ifndef HALF_LJ
 688     wco_vdw_S2  = gmx_cmplt_pr(rsq_S2, rcvdw2_S);
 689     wco_vdw_S3  = gmx_cmplt_pr(rsq_S3, rcvdw2_S);
 690 #endif
 691 #else
 692     /* Same cut-off for Coulomb and VdW, reuse the registers */
 693 #define     wco_vdw_S0    wco_S0
 694 #define     wco_vdw_S1    wco_S1
 695 #define     wco_vdw_S2    wco_S2
 696 #define     wco_vdw_S3    wco_S3
 697 #endif
 698
 699 #ifndef LJ_COMB_LB
 700     rinvsix_S0  = gmx_mul_pr(rinvsq_S0, gmx_mul_pr(rinvsq_S0, rinvsq_S0));
 701     rinvsix_S1  = gmx_mul_pr(rinvsq_S1, gmx_mul_pr(rinvsq_S1, rinvsq_S1));
 702 #ifdef EXCL_FORCES
 703     rinvsix_S0  = gmx_blendzero_pr(rinvsix_S0, int_S0);
 704     rinvsix_S1  = gmx_blendzero_pr(rinvsix_S1, int_S1);
 705 #endif
 706 #ifndef HALF_LJ
 707     rinvsix_S2  = gmx_mul_pr(rinvsq_S2, gmx_mul_pr(rinvsq_S2, rinvsq_S2));
 708     rinvsix_S3  = gmx_mul_pr(rinvsq_S3, gmx_mul_pr(rinvsq_S3, rinvsq_S3));
 709 #ifdef EXCL_FORCES
 710     rinvsix_S2  = gmx_blendzero_pr(rinvsix_S2, int_S2);
 711     rinvsix_S3  = gmx_blendzero_pr(rinvsix_S3, int_S3);
 712 #endif
 713 #endif
 714 #ifdef VDW_CUTOFF_CHECK
 715     rinvsix_S0  = gmx_blendzero_pr(rinvsix_S0, wco_vdw_S0);
 716     rinvsix_S1  = gmx_blendzero_pr(rinvsix_S1, wco_vdw_S1);
 717 #ifndef HALF_LJ
 718     rinvsix_S2  = gmx_blendzero_pr(rinvsix_S2, wco_vdw_S2);
 719     rinvsix_S3  = gmx_blendzero_pr(rinvsix_S3, wco_vdw_S3);
 720 #endif
 721 #endif
 722     FrLJ6_S0    = gmx_mul_pr(c6_S0, rinvsix_S0);
 723     FrLJ6_S1    = gmx_mul_pr(c6_S1, rinvsix_S1);
 724 #ifndef HALF_LJ
 725     FrLJ6_S2    = gmx_mul_pr(c6_S2, rinvsix_S2);
 726     FrLJ6_S3    = gmx_mul_pr(c6_S3, rinvsix_S3);
 727 #endif
 728     FrLJ12_S0   = gmx_mul_pr(c12_S0, gmx_mul_pr(rinvsix_S0, rinvsix_S0));
 729     FrLJ12_S1   = gmx_mul_pr(c12_S1, gmx_mul_pr(rinvsix_S1, rinvsix_S1));
 730 #ifndef HALF_LJ
 731     FrLJ12_S2   = gmx_mul_pr(c12_S2, gmx_mul_pr(rinvsix_S2, rinvsix_S2));
 732     FrLJ12_S3   = gmx_mul_pr(c12_S3, gmx_mul_pr(rinvsix_S3, rinvsix_S3));
 733 #endif
 734 #endif /* not LJ_COMB_LB */
 735
 736 #ifdef LJ_COMB_LB
 737     sir_S0      = gmx_mul_pr(sig_S0, rinv_S0);
 738     sir_S1      = gmx_mul_pr(sig_S1, rinv_S1);
 739 #ifndef HALF_LJ
 740     sir_S2      = gmx_mul_pr(sig_S2, rinv_S2);
 741     sir_S3      = gmx_mul_pr(sig_S3, rinv_S3);
 742 #endif
 743     sir2_S0     = gmx_mul_pr(sir_S0, sir_S0);
 744     sir2_S1     = gmx_mul_pr(sir_S1, sir_S1);
 745 #ifndef HALF_LJ
 746     sir2_S2     = gmx_mul_pr(sir_S2, sir_S2);
 747     sir2_S3     = gmx_mul_pr(sir_S3, sir_S3);
 748 #endif
 749     sir6_S0     = gmx_mul_pr(sir2_S0, gmx_mul_pr(sir2_S0, sir2_S0));
 750     sir6_S1     = gmx_mul_pr(sir2_S1, gmx_mul_pr(sir2_S1, sir2_S1));
 751 #ifdef EXCL_FORCES
 752     sir6_S0     = gmx_blendzero_pr(sir6_S0, int_S0);
 753     sir6_S1     = gmx_blendzero_pr(sir6_S1, int_S1);
 754 #endif
 755 #ifndef HALF_LJ
 756     sir6_S2     = gmx_mul_pr(sir2_S2, gmx_mul_pr(sir2_S2, sir2_S2));
 757     sir6_S3     = gmx_mul_pr(sir2_S3, gmx_mul_pr(sir2_S3, sir2_S3));
 758 #ifdef EXCL_FORCES
 759     sir6_S2     = gmx_blendzero_pr(sir6_S2, int_S2);
 760     sir6_S3     = gmx_blendzero_pr(sir6_S3, int_S3);
 761 #endif
 762 #endif
 763 #ifdef VDW_CUTOFF_CHECK
 764     sir6_S0     = gmx_blendzero_pr(sir6_S0, wco_vdw_S0);
 765     sir6_S1     = gmx_blendzero_pr(sir6_S1, wco_vdw_S1);
 766 #ifndef HALF_LJ
 767     sir6_S2     = gmx_blendzero_pr(sir6_S2, wco_vdw_S2);
 768     sir6_S3     = gmx_blendzero_pr(sir6_S3, wco_vdw_S3);
 769 #endif
 770 #endif
 771     FrLJ6_S0    = gmx_mul_pr(eps_S0, sir6_S0);
 772     FrLJ6_S1    = gmx_mul_pr(eps_S1, sir6_S1);
 773 #ifndef HALF_LJ
 774     FrLJ6_S2    = gmx_mul_pr(eps_S2, sir6_S2);
 775     FrLJ6_S3    = gmx_mul_pr(eps_S3, sir6_S3);
 776 #endif
 777     FrLJ12_S0   = gmx_mul_pr(FrLJ6_S0, sir6_S0);
 778     FrLJ12_S1   = gmx_mul_pr(FrLJ6_S1, sir6_S1);
 779 #ifndef HALF_LJ
 780     FrLJ12_S2   = gmx_mul_pr(FrLJ6_S2, sir6_S2);
 781     FrLJ12_S3   = gmx_mul_pr(FrLJ6_S3, sir6_S3);
 782 #endif
 783 #if defined CALC_ENERGIES
 784     /* We need C6 and C12 to calculate the LJ potential shift */
 785     sig2_S0     = gmx_mul_pr(sig_S0, sig_S0);
 786     sig2_S1     = gmx_mul_pr(sig_S1, sig_S1);
 787 #ifndef HALF_LJ
 788     sig2_S2     = gmx_mul_pr(sig_S2, sig_S2);
 789     sig2_S3     = gmx_mul_pr(sig_S3, sig_S3);
 790 #endif
 791     sig6_S0     = gmx_mul_pr(sig2_S0, gmx_mul_pr(sig2_S0, sig2_S0));
 792     sig6_S1     = gmx_mul_pr(sig2_S1, gmx_mul_pr(sig2_S1, sig2_S1));
 793 #ifndef HALF_LJ
 794     sig6_S2     = gmx_mul_pr(sig2_S2, gmx_mul_pr(sig2_S2, sig2_S2));
 795     sig6_S3     = gmx_mul_pr(sig2_S3, gmx_mul_pr(sig2_S3, sig2_S3));
 796 #endif
 797     c6_S0       = gmx_mul_pr(eps_S0, sig6_S0);
 798     c6_S1       = gmx_mul_pr(eps_S1, sig6_S1);
 799 #ifndef HALF_LJ
 800     c6_S2       = gmx_mul_pr(eps_S2, sig6_S2);
 801     c6_S3       = gmx_mul_pr(eps_S3, sig6_S3);
 802 #endif
 803     c12_S0      = gmx_mul_pr(c6_S0, sig6_S0);
 804     c12_S1      = gmx_mul_pr(c6_S1, sig6_S1);
 805 #ifndef HALF_LJ
 806     c12_S2      = gmx_mul_pr(c6_S2, sig6_S2);
 807     c12_S3      = gmx_mul_pr(c6_S3, sig6_S3);
 808 #endif
 809 #endif
 810 #endif /* LJ_COMB_LB */
 811
 812 #endif /* CALC_LJ */
 813
 814 #ifdef CALC_ENERGIES
 815 #ifdef ENERGY_GROUPS
 816     /* Extract the group pair index per j pair.
 817      * Energy groups are stored per i-cluster, so things get
 818      * complicated when the i- and j-cluster size don't match.
 819      */
 820     {
 821         int egps_j;
 822 #if UNROLLJ == 2
 823         egps_j    = nbat->energrp[cj>>1];
 824         egp_jj[0] = ((egps_j >> ((cj & 1)*egps_jshift)) & egps_jmask)*egps_jstride;
 825 #else
 826         /* We assume UNROLLI <= UNROLLJ */
 827         int jdi;
 828         for (jdi = 0; jdi < UNROLLJ/UNROLLI; jdi++)
 829         {
 830             int jj;
 831             egps_j = nbat->energrp[cj*(UNROLLJ/UNROLLI)+jdi];
 832             for (jj = 0; jj < (UNROLLI/2); jj++)
 833             {
 834                 egp_jj[jdi*(UNROLLI/2)+jj] = ((egps_j >> (jj*egps_jshift)) & egps_jmask)*egps_jstride;
 835             }
 836         }
 837 #endif
 838     }
 839 #endif
 840
 841 #ifdef CALC_COULOMB
 842 #ifndef ENERGY_GROUPS
 843     vctot_S      = gmx_add_pr(vctot_S, gmx_sum4_pr(vcoul_S0, vcoul_S1, vcoul_S2, vcoul_S3));
 844 #else
 845     add_ener_grp(vcoul_S0, vctp[0], egp_jj);
 846     add_ener_grp(vcoul_S1, vctp[1], egp_jj);
 847     add_ener_grp(vcoul_S2, vctp[2], egp_jj);
 848     add_ener_grp(vcoul_S3, vctp[3], egp_jj);
 849 #endif
 850 #endif
 851
 852 #ifdef CALC_LJ
 853     /* Calculate the LJ energies */
 854     VLJ6_S0     = gmx_mul_pr(sixth_S, gmx_sub_pr(FrLJ6_S0, gmx_mul_pr(c6_S0, sh_invrc6_S)));
 855     VLJ6_S1     = gmx_mul_pr(sixth_S, gmx_sub_pr(FrLJ6_S1, gmx_mul_pr(c6_S1, sh_invrc6_S)));
 856 #ifndef HALF_LJ
 857     VLJ6_S2     = gmx_mul_pr(sixth_S, gmx_sub_pr(FrLJ6_S2, gmx_mul_pr(c6_S2, sh_invrc6_S)));
 858     VLJ6_S3     = gmx_mul_pr(sixth_S, gmx_sub_pr(FrLJ6_S3, gmx_mul_pr(c6_S3, sh_invrc6_S)));
 859 #endif
 860     VLJ12_S0    = gmx_mul_pr(twelveth_S, gmx_sub_pr(FrLJ12_S0, gmx_mul_pr(c12_S0, sh_invrc12_S)));
 861     VLJ12_S1    = gmx_mul_pr(twelveth_S, gmx_sub_pr(FrLJ12_S1, gmx_mul_pr(c12_S1, sh_invrc12_S)));
 862 #ifndef HALF_LJ
 863     VLJ12_S2    = gmx_mul_pr(twelveth_S, gmx_sub_pr(FrLJ12_S2, gmx_mul_pr(c12_S2, sh_invrc12_S)));
 864     VLJ12_S3    = gmx_mul_pr(twelveth_S, gmx_sub_pr(FrLJ12_S3, gmx_mul_pr(c12_S3, sh_invrc12_S)));
 865 #endif
 866
 867     VLJ_S0      = gmx_sub_pr(VLJ12_S0, VLJ6_S0);
 868     VLJ_S1      = gmx_sub_pr(VLJ12_S1, VLJ6_S1);
 869 #ifndef HALF_LJ
 870     VLJ_S2      = gmx_sub_pr(VLJ12_S2, VLJ6_S2);
 871     VLJ_S3      = gmx_sub_pr(VLJ12_S3, VLJ6_S3);
 872 #endif
 873     /* The potential shift should be removed for pairs beyond cut-off */
 874     VLJ_S0      = gmx_blendzero_pr(VLJ_S0, wco_vdw_S0);
 875     VLJ_S1      = gmx_blendzero_pr(VLJ_S1, wco_vdw_S1);
 876 #ifndef HALF_LJ
 877     VLJ_S2      = gmx_blendzero_pr(VLJ_S2, wco_vdw_S2);
 878     VLJ_S3      = gmx_blendzero_pr(VLJ_S3, wco_vdw_S3);
 879 #endif
 880 #ifdef CHECK_EXCLS
 881     /* The potential shift should be removed for excluded pairs */
 882     VLJ_S0      = gmx_blendzero_pr(VLJ_S0, int_S0);
 883     VLJ_S1      = gmx_blendzero_pr(VLJ_S1, int_S1);
 884 #ifndef HALF_LJ
 885     VLJ_S2      = gmx_blendzero_pr(VLJ_S2, int_S2);
 886     VLJ_S3      = gmx_blendzero_pr(VLJ_S3, int_S3);
 887 #endif
 888 #endif
 889 #ifndef ENERGY_GROUPS
 890     Vvdwtot_S   = gmx_add_pr(Vvdwtot_S,
 891 #ifndef HALF_LJ
 892                              gmx_sum4_pr(VLJ_S0, VLJ_S1, VLJ_S2, VLJ_S3)
 893 #else
 894                              gmx_add_pr(VLJ_S0, VLJ_S1)
 895 #endif
 896                              );
 897 #else
 898     add_ener_grp(VLJ_S0, vvdwtp[0], egp_jj);
 899     add_ener_grp(VLJ_S1, vvdwtp[1], egp_jj);
 900 #ifndef HALF_LJ
 901     add_ener_grp(VLJ_S2, vvdwtp[2], egp_jj);
 902     add_ener_grp(VLJ_S3, vvdwtp[3], egp_jj);
 903 #endif
 904 #endif
 905 #endif /* CALC_LJ */
 906 #endif /* CALC_ENERGIES */
 907
 908 #ifdef CALC_LJ
 909     fscal_S0    = gmx_mul_pr(rinvsq_S0,
 910 #ifdef CALC_COULOMB
 911                                gmx_add_pr(frcoul_S0,
 912 #else
 913                                (
 914 #endif
 915                                           gmx_sub_pr(FrLJ12_S0, FrLJ6_S0)));
 916     fscal_S1    = gmx_mul_pr(rinvsq_S1,
 917 #ifdef CALC_COULOMB
 918                                gmx_add_pr(frcoul_S1,
 919 #else
 920                                (
 921 #endif
 922                                           gmx_sub_pr(FrLJ12_S1, FrLJ6_S1)));
 923 #else
 924     fscal_S0    = gmx_mul_pr(rinvsq_S0, frcoul_S0);
 925     fscal_S1    = gmx_mul_pr(rinvsq_S1, frcoul_S1);
 926 #endif /* CALC_LJ */
 927 #if defined CALC_LJ && !defined HALF_LJ
 928     fscal_S2    = gmx_mul_pr(rinvsq_S2,
 929 #ifdef CALC_COULOMB
 930                                gmx_add_pr(frcoul_S2,
 931 #else
 932                                (
 933 #endif
 934                                           gmx_sub_pr(FrLJ12_S2, FrLJ6_S2)));
 935     fscal_S3    = gmx_mul_pr(rinvsq_S3,
 936 #ifdef CALC_COULOMB
 937                                gmx_add_pr(frcoul_S3,
 938 #else
 939                                (
 940 #endif
 941                                           gmx_sub_pr(FrLJ12_S3, FrLJ6_S3)));
 942 #else
 943     /* Atom 2 and 3 don't have LJ, so only add Coulomb forces */
 944     fscal_S2    = gmx_mul_pr(rinvsq_S2, frcoul_S2);
 945     fscal_S3    = gmx_mul_pr(rinvsq_S3, frcoul_S3);
 946 #endif
 947
 948     /* Calculate temporary vectorial force */
 949     tx_S0       = gmx_mul_pr(fscal_S0, dx_S0);
 950     tx_S1       = gmx_mul_pr(fscal_S1, dx_S1);
 951     tx_S2       = gmx_mul_pr(fscal_S2, dx_S2);
 952     tx_S3       = gmx_mul_pr(fscal_S3, dx_S3);
 953     ty_S0       = gmx_mul_pr(fscal_S0, dy_S0);
 954     ty_S1       = gmx_mul_pr(fscal_S1, dy_S1);
 955     ty_S2       = gmx_mul_pr(fscal_S2, dy_S2);
 956     ty_S3       = gmx_mul_pr(fscal_S3, dy_S3);
 957     tz_S0       = gmx_mul_pr(fscal_S0, dz_S0);
 958     tz_S1       = gmx_mul_pr(fscal_S1, dz_S1);
 959     tz_S2       = gmx_mul_pr(fscal_S2, dz_S2);
 960     tz_S3       = gmx_mul_pr(fscal_S3, dz_S3);
 961
 962     /* Increment i atom force */
 963     fix_S0      = gmx_add_pr(fix_S0, tx_S0);
 964     fix_S1      = gmx_add_pr(fix_S1, tx_S1);
 965     fix_S2      = gmx_add_pr(fix_S2, tx_S2);
 966     fix_S3      = gmx_add_pr(fix_S3, tx_S3);
 967     fiy_S0      = gmx_add_pr(fiy_S0, ty_S0);
 968     fiy_S1      = gmx_add_pr(fiy_S1, ty_S1);
 969     fiy_S2      = gmx_add_pr(fiy_S2, ty_S2);
 970     fiy_S3      = gmx_add_pr(fiy_S3, ty_S3);
 971     fiz_S0      = gmx_add_pr(fiz_S0, tz_S0);
 972     fiz_S1      = gmx_add_pr(fiz_S1, tz_S1);
 973     fiz_S2      = gmx_add_pr(fiz_S2, tz_S2);
 974     fiz_S3      = gmx_add_pr(fiz_S3, tz_S3);
 975
 976     /* Decrement j atom force */
 977     gmx_store_pr(f+ajx,
 978                  gmx_sub_pr( gmx_load_pr(f+ajx), gmx_sum4_pr(tx_S0, tx_S1, tx_S2, tx_S3) ));
 979     gmx_store_pr(f+ajy,
 980                  gmx_sub_pr( gmx_load_pr(f+ajy), gmx_sum4_pr(ty_S0, ty_S1, ty_S2, ty_S3) ));
 981     gmx_store_pr(f+ajz,
 982                  gmx_sub_pr( gmx_load_pr(f+ajz), gmx_sum4_pr(tz_S0, tz_S1, tz_S2, tz_S3) ));
 983 }
 984
 985 #undef  rinv_ex_S0
 986 #undef  rinv_ex_S1
 987 #undef  rinv_ex_S2
 988 #undef  rinv_ex_S3
 989
 990 #undef  wco_vdw_S0
 991 #undef  wco_vdw_S1
 992 #undef  wco_vdw_S2
 993 #undef  wco_vdw_S3
 994
 995 #undef  NBNXN_CUTOFF_USE_BLENDV
 996
 997 #undef  EXCL_FORCES