src/gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_inner.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 /* This is the innermost loop contents for the 4 x N atom SIMD kernel.
  37  * This flavor of the kernel duplicates the data for N j-particles in
  38  * 2xN wide SIMD registers to do operate on 2 i-particles at once.
  39  * This leads to 4/2=2 sets of most instructions. Therefore we call
  40  * this kernel 2x(N+N) = 2xnn
  41  *
  42  * This 2xnn kernel is basically the 4xn equivalent with half the registers
  43  * and instructions removed.
  44  *
  45  * An alternative would be to load to different cluster of N j-particles
  46  * into SIMD registers, giving a 4x(N+N) kernel. This doubles the amount
  47  * of instructions, which could lead to better scheduling. But we actually
  48  * observed worse scheduling for the AVX-256 4x8 normal analytical PME
  49  * kernel, which has a lower pair throughput than 2x(4+4) with gcc 4.7.
  50  * It could be worth trying this option, but it takes some more effort.
  51  * This 2xnn kernel is basically the 4xn equivalent with
  52  */
  53
  54
  55 /* When calculating RF or Ewald interactions we calculate the electrostatic/LJ
  56  * forces on excluded atom pairs here in the non-bonded loops.
  57  * But when energies and/or virial is required we calculate them
  58  * separately to as then it is easier to separate the energy and virial
  59  * contributions.
  60  */
  61 #if defined CHECK_EXCLS && (defined CALC_COULOMB || defined LJ_EWALD_GEOM)
  62 #define EXCL_FORCES
  63 #endif
  64
  65 /* Without exclusions and energies we only need to mask the cut-off,
  66  * this can be faster with blendv.
  67  */
  68 #if !(defined CHECK_EXCLS || defined CALC_ENERGIES) && defined GMX_SIMD_HAVE_BLENDV && !defined COUNT_PAIRS
  69 /* With RF and tabulated Coulomb we replace cmp+and with sub+blendv.
  70  * With gcc this is slower, except for RF on Sandy Bridge.
  71  * Tested with gcc 4.6.2, 4.6.3 and 4.7.1.
  72  */
  73 #if (defined CALC_COUL_RF || defined CALC_COUL_TAB) && (!defined __GNUC__ || (defined CALC_COUL_RF && defined GMX_SIMD_X86_AVX_256_OR_HIGHER))
  74 #define NBNXN_CUTOFF_USE_BLENDV
  75 #endif
  76 /* With analytical Ewald we replace cmp+and+and with sub+blendv+blendv.
  77  * This is only faster with icc on Sandy Bridge (PS kernel slower than gcc 4.7).
  78  * Tested with icc 13.
  79  */
  80 #if defined CALC_COUL_EWALD && defined __INTEL_COMPILER && defined GMX_SIMD_X86_AVX_256_OR_HIGHER
  81 #define NBNXN_CUTOFF_USE_BLENDV
  82 #endif
  83 #endif
  84
  85 {
  86     int        cj, aj, ajx, ajy, ajz;
  87
  88 #ifdef ENERGY_GROUPS
  89     /* Energy group indices for two atoms packed into one int */
  90     int        egp_jj[UNROLLJ/2];
  91 #endif
  92
  93 #ifdef CHECK_EXCLS
  94     /* Interaction (non-exclusion) mask of all 1's or 0's */
  95     gmx_simd_bool_t  interact_S0;
  96     gmx_simd_bool_t  interact_S2;
  97 #endif
  98
  99     gmx_simd_real_t  jx_S, jy_S, jz_S;
 100     gmx_simd_real_t  dx_S0, dy_S0, dz_S0;
 101     gmx_simd_real_t  dx_S2, dy_S2, dz_S2;
 102     gmx_simd_real_t  tx_S0, ty_S0, tz_S0;
 103     gmx_simd_real_t  tx_S2, ty_S2, tz_S2;
 104     gmx_simd_real_t  rsq_S0, rinv_S0, rinvsq_S0;
 105     gmx_simd_real_t  rsq_S2, rinv_S2, rinvsq_S2;
 106 #ifndef NBNXN_CUTOFF_USE_BLENDV
 107     /* wco: within cut-off, mask of all 1's or 0's */
 108     gmx_simd_bool_t  wco_S0;
 109     gmx_simd_bool_t  wco_S2;
 110 #endif
 111 #ifdef VDW_CUTOFF_CHECK
 112     gmx_simd_bool_t  wco_vdw_S0;
 113 #ifndef HALF_LJ
 114     gmx_simd_bool_t  wco_vdw_S2;
 115 #endif
 116 #endif
 117
 118 #if (defined CALC_COULOMB && defined CALC_COUL_TAB) || defined LJ_FORCE_SWITCH || defined LJ_POT_SWITCH
 119     gmx_simd_real_t r_S0;
 120     gmx_simd_real_t r_S2;
 121 #endif
 122
 123 #if defined LJ_FORCE_SWITCH || defined LJ_POT_SWITCH
 124     gmx_simd_real_t  rsw_S0, rsw2_S0, rsw2_r_S0;
 125 #ifndef HALF_LJ
 126     gmx_simd_real_t  rsw_S2, rsw2_S2, rsw2_r_S2;
 127 #endif
 128 #endif
 129
 130 #ifdef CALC_COULOMB
 131 #ifdef CHECK_EXCLS
 132     /* 1/r masked with the interaction mask */
 133     gmx_simd_real_t  rinv_ex_S0;
 134     gmx_simd_real_t  rinv_ex_S2;
 135 #endif
 136     gmx_simd_real_t  jq_S;
 137     gmx_simd_real_t  qq_S0;
 138     gmx_simd_real_t  qq_S2;
 139 #ifdef CALC_COUL_TAB
 140     /* The force (PME mesh force) we need to subtract from 1/r^2 */
 141     gmx_simd_real_t  fsub_S0;
 142     gmx_simd_real_t  fsub_S2;
 143 #endif
 144 #ifdef CALC_COUL_EWALD
 145     gmx_simd_real_t  brsq_S0, brsq_S2;
 146     gmx_simd_real_t  ewcorr_S0, ewcorr_S2;
 147 #endif
 148
 149     /* frcoul = (1/r - fsub)*r */
 150     gmx_simd_real_t  frcoul_S0;
 151     gmx_simd_real_t  frcoul_S2;
 152 #ifdef CALC_COUL_TAB
 153     /* For tables: r, rs=r/sp, rf=floor(rs), frac=rs-rf */
 154     gmx_simd_real_t         rs_S0, rf_S0, frac_S0;
 155     gmx_simd_real_t         rs_S2, rf_S2, frac_S2;
 156     /* Table index: rs truncated to an int */
 157     gmx_simd_int32_t        ti_S0, ti_S2;
 158     /* Linear force table values */
 159     gmx_simd_real_t         ctab0_S0, ctab1_S0;
 160     gmx_simd_real_t         ctab0_S2, ctab1_S2;
 161 #ifdef CALC_ENERGIES
 162     /* Quadratic energy table value */
 163     gmx_simd_real_t  ctabv_S0;
 164     gmx_simd_real_t  ctabv_S2;
 165 #endif
 166 #endif
 167 #if defined CALC_ENERGIES && (defined CALC_COUL_EWALD || defined CALC_COUL_TAB)
 168     /* The potential (PME mesh) we need to subtract from 1/r */
 169     gmx_simd_real_t  vc_sub_S0;
 170     gmx_simd_real_t  vc_sub_S2;
 171 #endif
 172 #ifdef CALC_ENERGIES
 173     /* Electrostatic potential */
 174     gmx_simd_real_t  vcoul_S0;
 175     gmx_simd_real_t  vcoul_S2;
 176 #endif
 177 #endif
 178     /* The force times 1/r */
 179     gmx_simd_real_t  fscal_S0;
 180     gmx_simd_real_t  fscal_S2;
 181
 182 #ifdef CALC_LJ
 183 #ifdef LJ_COMB_LB
 184     /* LJ sigma_j/2 and sqrt(epsilon_j) */
 185     gmx_simd_real_t  hsig_j_S, seps_j_S;
 186     /* LJ sigma_ij and epsilon_ij */
 187     gmx_simd_real_t  sig_S0, eps_S0;
 188 #ifndef HALF_LJ
 189     gmx_simd_real_t  sig_S2, eps_S2;
 190 #endif
 191 #ifdef CALC_ENERGIES
 192     gmx_simd_real_t  sig2_S0, sig6_S0;
 193 #ifndef HALF_LJ
 194     gmx_simd_real_t  sig2_S2, sig6_S2;
 195 #endif
 196 #endif /* LJ_COMB_LB */
 197 #endif /* CALC_LJ */
 198
 199 #ifdef LJ_COMB_GEOM
 200     gmx_simd_real_t  c6s_j_S, c12s_j_S;
 201 #endif
 202
 203 #if defined LJ_COMB_GEOM || defined LJ_COMB_LB || defined LJ_EWALD_GEOM
 204     /* Index for loading LJ parameters, complicated when interleaving */
 205     int         aj2;
 206 #endif
 207
 208 #ifndef FIX_LJ_C
 209     /* LJ C6 and C12 parameters, used with geometric comb. rule */
 210     gmx_simd_real_t  c6_S0, c12_S0;
 211 #ifndef HALF_LJ
 212     gmx_simd_real_t  c6_S2, c12_S2;
 213 #endif
 214 #endif
 215
 216     /* Intermediate variables for LJ calculation */
 217 #ifndef LJ_COMB_LB
 218     gmx_simd_real_t  rinvsix_S0;
 219 #ifndef HALF_LJ
 220     gmx_simd_real_t  rinvsix_S2;
 221 #endif
 222 #endif
 223 #ifdef LJ_COMB_LB
 224     gmx_simd_real_t  sir_S0, sir2_S0, sir6_S0;
 225 #ifndef HALF_LJ
 226     gmx_simd_real_t  sir_S2, sir2_S2, sir6_S2;
 227 #endif
 228 #endif
 229
 230     gmx_simd_real_t  FrLJ6_S0, FrLJ12_S0, frLJ_S0;
 231 #ifndef HALF_LJ
 232     gmx_simd_real_t  FrLJ6_S2, FrLJ12_S2, frLJ_S2;
 233 #endif
 234 #if defined CALC_ENERGIES || defined LJ_POT_SWITCH
 235     gmx_simd_real_t  VLJ6_S0, VLJ12_S0, VLJ_S0;
 236 #ifndef HALF_LJ
 237     gmx_simd_real_t  VLJ6_S2, VLJ12_S2, VLJ_S2;
 238 #endif
 239 #endif
 240 #endif /* CALC_LJ */
 241
 242     gmx_mm_hpr fjx_S, fjy_S, fjz_S;
 243
 244     /* j-cluster index */
 245     cj            = l_cj[cjind].cj;
 246
 247     /* Atom indices (of the first atom in the cluster) */
 248     aj            = cj*UNROLLJ;
 249 #if defined CALC_LJ && (defined LJ_COMB_GEOM || defined LJ_COMB_LB || defined LJ_EWALD_GEOM)
 250 #if UNROLLJ == STRIDE
 251     aj2           = aj*2;
 252 #else
 253     aj2           = (cj>>1)*2*STRIDE + (cj & 1)*UNROLLJ;
 254 #endif
 255 #endif
 256 #if UNROLLJ == STRIDE
 257     ajx           = aj*DIM;
 258 #else
 259     ajx           = (cj>>1)*DIM*STRIDE + (cj & 1)*UNROLLJ;
 260 #endif
 261     ajy           = ajx + STRIDE;
 262     ajz           = ajy + STRIDE;
 263
 264 #ifdef CHECK_EXCLS
 265     gmx_load_simd_2xnn_interactions(l_cj[cjind].excl,
 266                                     filter_S0, filter_S2,
 267                                     &interact_S0, &interact_S2);
 268 #endif /* CHECK_EXCLS */
 269
 270     /* load j atom coordinates */
 271     gmx_loaddh_pr(&jx_S, x+ajx);
 272     gmx_loaddh_pr(&jy_S, x+ajy);
 273     gmx_loaddh_pr(&jz_S, x+ajz);
 274
 275     /* Calculate distance */
 276     dx_S0       = gmx_simd_sub_r(ix_S0, jx_S);
 277     dy_S0       = gmx_simd_sub_r(iy_S0, jy_S);
 278     dz_S0       = gmx_simd_sub_r(iz_S0, jz_S);
 279     dx_S2       = gmx_simd_sub_r(ix_S2, jx_S);
 280     dy_S2       = gmx_simd_sub_r(iy_S2, jy_S);
 281     dz_S2       = gmx_simd_sub_r(iz_S2, jz_S);
 282
 283     /* rsq = dx*dx+dy*dy+dz*dz */
 284     rsq_S0      = gmx_simd_calc_rsq_r(dx_S0, dy_S0, dz_S0);
 285     rsq_S2      = gmx_simd_calc_rsq_r(dx_S2, dy_S2, dz_S2);
 286
 287 #ifndef NBNXN_CUTOFF_USE_BLENDV
 288     wco_S0      = gmx_simd_cmplt_r(rsq_S0, rc2_S);
 289     wco_S2      = gmx_simd_cmplt_r(rsq_S2, rc2_S);
 290 #endif
 291
 292 #ifdef CHECK_EXCLS
 293 #ifdef EXCL_FORCES
 294     /* Only remove the (sub-)diagonal to avoid double counting */
 295 #if UNROLLJ == UNROLLI
 296     if (cj == ci_sh)
 297     {
 298         wco_S0  = gmx_simd_and_b(wco_S0, diagonal_mask_S0);
 299         wco_S2  = gmx_simd_and_b(wco_S2, diagonal_mask_S2);
 300     }
 301 #else
 302 #if UNROLLJ == 2*UNROLLI
 303     if (cj*2 == ci_sh)
 304     {
 305         wco_S0  = gmx_simd_and_b(wco_S0, diagonal_mask0_S0);
 306         wco_S2  = gmx_simd_and_b(wco_S2, diagonal_mask0_S2);
 307     }
 308     else if (cj*2 + 1 == ci_sh)
 309     {
 310         wco_S0  = gmx_simd_and_b(wco_S0, diagonal_mask1_S0);
 311         wco_S2  = gmx_simd_and_b(wco_S2, diagonal_mask1_S2);
 312     }
 313 #else
 314 #error "only UNROLLJ == UNROLLI*(1 or 2) currently supported in 2xnn kernels"
 315 #endif
 316 #endif
 317 #else /* EXCL_FORCES */
 318       /* No exclusion forces: remove all excluded atom pairs from the list */
 319     wco_S0      = gmx_simd_and_b(wco_S0, interact_S0);
 320     wco_S2      = gmx_simd_and_b(wco_S2, interact_S2);
 321 #endif
 322 #endif
 323
 324 #ifdef COUNT_PAIRS
 325     {
 326         int  i, j;
 327         real tmpa[2*GMX_SIMD_REAL_WIDTH], *tmp;
 328         tmp = gmx_simd_align_r(tmpa);
 329         for (i = 0; i < UNROLLI; i += 2)
 330         {
 331             gmx_simd_store_r(tmp, i == 0 ? wco_S0 : wco_S2);
 332             for (j = 0; j < 2*UNROLLJ; j++)
 333             {
 334                 if (!(tmp[j] == 0))
 335                 {
 336                     npair++;
 337                 }
 338             }
 339         }
 340     }
 341 #endif
 342
 343 #ifdef CHECK_EXCLS
 344     /* For excluded pairs add a small number to avoid r^-6 = NaN */
 345     rsq_S0      = gmx_simd_add_r(rsq_S0, gmx_simd_blendnotzero_r(avoid_sing_S, interact_S0));
 346     rsq_S2      = gmx_simd_add_r(rsq_S2, gmx_simd_blendnotzero_r(avoid_sing_S, interact_S2));
 347 #endif
 348
 349     /* Calculate 1/r */
 350     rinv_S0     = gmx_simd_invsqrt_r(rsq_S0);
 351     rinv_S2     = gmx_simd_invsqrt_r(rsq_S2);
 352
 353 #ifdef CALC_COULOMB
 354     /* Load parameters for j atom */
 355     gmx_loaddh_pr(&jq_S, q+aj);
 356     qq_S0       = gmx_simd_mul_r(iq_S0, jq_S);
 357     qq_S2       = gmx_simd_mul_r(iq_S2, jq_S);
 358 #endif
 359
 360 #ifdef CALC_LJ
 361
 362 #if !defined LJ_COMB_GEOM && !defined LJ_COMB_LB && !defined FIX_LJ_C
 363     load_lj_pair_params2(nbfp0, nbfp1, type, aj, &c6_S0, &c12_S0);
 364 #ifndef HALF_LJ
 365     load_lj_pair_params2(nbfp2, nbfp3, type, aj, &c6_S2, &c12_S2);
 366 #endif
 367 #endif /* not defined any LJ rule */
 368
 369 #ifdef LJ_COMB_GEOM
 370     gmx_loaddh_pr(&c6s_j_S,  ljc+aj2+0);
 371     gmx_loaddh_pr(&c12s_j_S, ljc+aj2+STRIDE);
 372     c6_S0       = gmx_simd_mul_r(c6s_S0, c6s_j_S );
 373 #ifndef HALF_LJ
 374     c6_S2       = gmx_simd_mul_r(c6s_S2, c6s_j_S );
 375 #endif
 376     c12_S0      = gmx_simd_mul_r(c12s_S0, c12s_j_S);
 377 #ifndef HALF_LJ
 378     c12_S2      = gmx_simd_mul_r(c12s_S2, c12s_j_S);
 379 #endif
 380 #endif /* LJ_COMB_GEOM */
 381
 382 #ifdef LJ_COMB_LB
 383     gmx_loaddh_pr(&hsig_j_S, ljc+aj2+0);
 384     gmx_loaddh_pr(&seps_j_S, ljc+aj2+STRIDE);
 385
 386     sig_S0      = gmx_simd_add_r(hsig_i_S0, hsig_j_S);
 387     eps_S0      = gmx_simd_mul_r(seps_i_S0, seps_j_S);
 388 #ifndef HALF_LJ
 389     sig_S2      = gmx_simd_add_r(hsig_i_S2, hsig_j_S);
 390     eps_S2      = gmx_simd_mul_r(seps_i_S2, seps_j_S);
 391 #endif
 392 #endif /* LJ_COMB_LB */
 393
 394 #endif /* CALC_LJ */
 395
 396 #ifndef NBNXN_CUTOFF_USE_BLENDV
 397     rinv_S0     = gmx_simd_blendzero_r(rinv_S0, wco_S0);
 398     rinv_S2     = gmx_simd_blendzero_r(rinv_S2, wco_S2);
 399 #else
 400     /* This needs to be modified: It makes assumptions about the internal storage
 401      * of the SIMD representation, in particular that the blendv instruction always
 402      * selects based on the sign bit. If the performance is really critical, it
 403      * should be turned into a function that is platform-specific.
 404      */
 405     /* We only need to mask for the cut-off: blendv is faster */
 406     rinv_S0     = gmx_simd_blendv_r(rinv_S0, zero_S, gmx_simd_sub_r(rc2_S, rsq_S0));
 407     rinv_S2     = gmx_simd_blendv_r(rinv_S2, zero_S, gmx_simd_sub_r(rc2_S, rsq_S2));
 408 #endif
 409
 410     rinvsq_S0   = gmx_simd_mul_r(rinv_S0, rinv_S0);
 411     rinvsq_S2   = gmx_simd_mul_r(rinv_S2, rinv_S2);
 412
 413 #ifdef CALC_COULOMB
 414     /* Note that here we calculate force*r, not the usual force/r.
 415      * This allows avoiding masking the reaction-field contribution,
 416      * as frcoul is later multiplied by rinvsq which has been
 417      * masked with the cut-off check.
 418      */
 419
 420 #ifdef EXCL_FORCES
 421     /* Only add 1/r for non-excluded atom pairs */
 422     rinv_ex_S0  = gmx_simd_blendzero_r(rinv_S0, interact_S0);
 423     rinv_ex_S2  = gmx_simd_blendzero_r(rinv_S2, interact_S2);
 424 #else
 425     /* No exclusion forces, we always need 1/r */
 426 #define     rinv_ex_S0    rinv_S0
 427 #define     rinv_ex_S2    rinv_S2
 428 #endif
 429
 430 #ifdef CALC_COUL_RF
 431     /* Electrostatic interactions */
 432     frcoul_S0   = gmx_simd_mul_r(qq_S0, gmx_simd_fmadd_r(rsq_S0, mrc_3_S, rinv_ex_S0));
 433     frcoul_S2   = gmx_simd_mul_r(qq_S2, gmx_simd_fmadd_r(rsq_S2, mrc_3_S, rinv_ex_S2));
 434
 435 #ifdef CALC_ENERGIES
 436     vcoul_S0    = gmx_simd_mul_r(qq_S0, gmx_simd_add_r(rinv_ex_S0, gmx_simd_add_r(gmx_simd_mul_r(rsq_S0, hrc_3_S), moh_rc_S)));
 437     vcoul_S2    = gmx_simd_mul_r(qq_S2, gmx_simd_add_r(rinv_ex_S2, gmx_simd_add_r(gmx_simd_mul_r(rsq_S2, hrc_3_S), moh_rc_S)));
 438 #endif
 439 #endif
 440
 441 #ifdef CALC_COUL_EWALD
 442     /* We need to mask (or limit) rsq for the cut-off,
 443      * as large distances can cause an overflow in gmx_pmecorrF/V.
 444      */
 445 #ifndef NBNXN_CUTOFF_USE_BLENDV
 446     brsq_S0     = gmx_simd_mul_r(beta2_S, gmx_simd_blendzero_r(rsq_S0, wco_S0));
 447     brsq_S2     = gmx_simd_mul_r(beta2_S, gmx_simd_blendzero_r(rsq_S2, wco_S2));
 448 #else
 449     /* Strangely, putting mul on a separate line is slower (icc 13) */
 450     brsq_S0     = gmx_simd_mul_r(beta2_S, gmx_simd_blendv_r(rsq_S0, zero_S, gmx_simd_sub_r(rc2_S, rsq_S0)));
 451     brsq_S2     = gmx_simd_mul_r(beta2_S, gmx_simd_blendv_r(rsq_S2, zero_S, gmx_simd_sub_r(rc2_S, rsq_S2)));
 452 #endif
 453     ewcorr_S0   = gmx_simd_mul_r(gmx_simd_pmecorrF_r(brsq_S0), beta_S);
 454     ewcorr_S2   = gmx_simd_mul_r(gmx_simd_pmecorrF_r(brsq_S2), beta_S);
 455     frcoul_S0   = gmx_simd_mul_r(qq_S0, gmx_simd_fmadd_r(ewcorr_S0, brsq_S0, rinv_ex_S0));
 456     frcoul_S2   = gmx_simd_mul_r(qq_S2, gmx_simd_fmadd_r(ewcorr_S2, brsq_S2, rinv_ex_S2));
 457
 458 #ifdef CALC_ENERGIES
 459     vc_sub_S0   = gmx_simd_mul_r(gmx_simd_pmecorrV_r(brsq_S0), beta_S);
 460     vc_sub_S2   = gmx_simd_mul_r(gmx_simd_pmecorrV_r(brsq_S2), beta_S);
 461 #endif
 462
 463 #endif /* CALC_COUL_EWALD */
 464
 465 #ifdef CALC_COUL_TAB
 466     /* Electrostatic interactions */
 467     r_S0        = gmx_simd_mul_r(rsq_S0, rinv_S0);
 468     r_S2        = gmx_simd_mul_r(rsq_S2, rinv_S2);
 469     /* Convert r to scaled table units */
 470     rs_S0       = gmx_simd_mul_r(r_S0, invtsp_S);
 471     rs_S2       = gmx_simd_mul_r(r_S2, invtsp_S);
 472     /* Truncate scaled r to an int */
 473     ti_S0       = gmx_simd_cvtt_r2i(rs_S0);
 474     ti_S2       = gmx_simd_cvtt_r2i(rs_S2);
 475 #ifdef GMX_SIMD_HAVE_TRUNC
 476     rf_S0       = gmx_simd_trunc_r(rs_S0);
 477     rf_S2       = gmx_simd_trunc_r(rs_S2);
 478 #else
 479     rf_S0       = gmx_simd_cvt_i2r(ti_S0);
 480     rf_S2       = gmx_simd_cvt_i2r(ti_S2);
 481 #endif
 482     frac_S0     = gmx_simd_sub_r(rs_S0, rf_S0);
 483     frac_S2     = gmx_simd_sub_r(rs_S2, rf_S2);
 484
 485     /* Load and interpolate table forces and possibly energies.
 486      * Force and energy can be combined in one table, stride 4: FDV0
 487      * or in two separate tables with stride 1: F and V
 488      * Currently single precision uses FDV0, double F and V.
 489      */
 490 #ifndef CALC_ENERGIES
 491     load_table_f(tab_coul_F, ti_S0, ti0, &ctab0_S0, &ctab1_S0);
 492     load_table_f(tab_coul_F, ti_S2, ti2, &ctab0_S2, &ctab1_S2);
 493 #else
 494 #ifdef TAB_FDV0
 495     load_table_f_v(tab_coul_F, ti_S0, ti0, &ctab0_S0, &ctab1_S0, &ctabv_S0);
 496     load_table_f_v(tab_coul_F, ti_S2, ti2, &ctab0_S2, &ctab1_S2, &ctabv_S2);
 497 #else
 498     load_table_f_v(tab_coul_F, tab_coul_V, ti_S0, ti0, &ctab0_S0, &ctab1_S0, &ctabv_S0);
 499     load_table_f_v(tab_coul_F, tab_coul_V, ti_S2, ti2, &ctab0_S2, &ctab1_S2, &ctabv_S2);
 500 #endif
 501 #endif
 502     fsub_S0     = gmx_simd_add_r(ctab0_S0, gmx_simd_mul_r(frac_S0, ctab1_S0));
 503     fsub_S2     = gmx_simd_add_r(ctab0_S2, gmx_simd_mul_r(frac_S2, ctab1_S2));
 504     frcoul_S0   = gmx_simd_mul_r(qq_S0, gmx_simd_sub_r(rinv_ex_S0, gmx_simd_mul_r(fsub_S0, r_S0)));
 505     frcoul_S2   = gmx_simd_mul_r(qq_S2, gmx_simd_sub_r(rinv_ex_S2, gmx_simd_mul_r(fsub_S2, r_S2)));
 506
 507 #ifdef CALC_ENERGIES
 508     vc_sub_S0   = gmx_simd_add_r(ctabv_S0, gmx_simd_mul_r(gmx_simd_mul_r(mhalfsp_S, frac_S0), gmx_simd_add_r(ctab0_S0, fsub_S0)));
 509     vc_sub_S2   = gmx_simd_add_r(ctabv_S2, gmx_simd_mul_r(gmx_simd_mul_r(mhalfsp_S, frac_S2), gmx_simd_add_r(ctab0_S2, fsub_S2)));
 510 #endif
 511 #endif /* CALC_COUL_TAB */
 512
 513 #if defined CALC_ENERGIES && (defined CALC_COUL_EWALD || defined CALC_COUL_TAB)
 514 #ifndef NO_SHIFT_EWALD
 515     /* Add Ewald potential shift to vc_sub for convenience */
 516 #ifdef CHECK_EXCLS
 517     vc_sub_S0   = gmx_simd_add_r(vc_sub_S0, gmx_simd_blendzero_r(sh_ewald_S, interact_S0));
 518     vc_sub_S2   = gmx_simd_add_r(vc_sub_S2, gmx_simd_blendzero_r(sh_ewald_S, interact_S2));
 519 #else
 520     vc_sub_S0   = gmx_simd_add_r(vc_sub_S0, sh_ewald_S);
 521     vc_sub_S2   = gmx_simd_add_r(vc_sub_S2, sh_ewald_S);
 522 #endif
 523 #endif
 524
 525     vcoul_S0    = gmx_simd_mul_r(qq_S0, gmx_simd_sub_r(rinv_ex_S0, vc_sub_S0));
 526     vcoul_S2    = gmx_simd_mul_r(qq_S2, gmx_simd_sub_r(rinv_ex_S2, vc_sub_S2));
 527 #endif
 528
 529 #ifdef CALC_ENERGIES
 530     /* Mask energy for cut-off and diagonal */
 531     vcoul_S0    = gmx_simd_blendzero_r(vcoul_S0, wco_S0);
 532     vcoul_S2    = gmx_simd_blendzero_r(vcoul_S2, wco_S2);
 533 #endif
 534
 535 #endif /* CALC_COULOMB */
 536
 537 #ifdef CALC_LJ
 538     /* Lennard-Jones interaction */
 539
 540 #ifdef VDW_CUTOFF_CHECK
 541     wco_vdw_S0  = gmx_simd_cmplt_r(rsq_S0, rcvdw2_S);
 542 #ifndef HALF_LJ
 543     wco_vdw_S2  = gmx_simd_cmplt_r(rsq_S2, rcvdw2_S);
 544 #endif
 545 #else
 546     /* Same cut-off for Coulomb and VdW, reuse the registers */
 547 #define     wco_vdw_S0    wco_S0
 548 #define     wco_vdw_S2    wco_S2
 549 #endif
 550
 551 #ifndef LJ_COMB_LB
 552     rinvsix_S0  = gmx_simd_mul_r(rinvsq_S0, gmx_simd_mul_r(rinvsq_S0, rinvsq_S0));
 553 #ifdef EXCL_FORCES
 554     rinvsix_S0  = gmx_simd_blendzero_r(rinvsix_S0, interact_S0);
 555 #endif
 556 #ifndef HALF_LJ
 557     rinvsix_S2  = gmx_simd_mul_r(rinvsq_S2, gmx_simd_mul_r(rinvsq_S2, rinvsq_S2));
 558 #ifdef EXCL_FORCES
 559     rinvsix_S2  = gmx_simd_blendzero_r(rinvsix_S2, interact_S2);
 560 #endif
 561 #endif
 562
 563 #if defined LJ_CUT || defined LJ_POT_SWITCH
 564     /* We have plain LJ or LJ-PME with simple C6/6 C12/12 coefficients */
 565     FrLJ6_S0    = gmx_simd_mul_r(c6_S0, rinvsix_S0);
 566 #ifndef HALF_LJ
 567     FrLJ6_S2    = gmx_simd_mul_r(c6_S2, rinvsix_S2);
 568 #endif
 569     FrLJ12_S0   = gmx_simd_mul_r(c12_S0, gmx_simd_mul_r(rinvsix_S0, rinvsix_S0));
 570 #ifndef HALF_LJ
 571     FrLJ12_S2   = gmx_simd_mul_r(c12_S2, gmx_simd_mul_r(rinvsix_S2, rinvsix_S2));
 572 #endif
 573 #endif
 574
 575 #if defined LJ_FORCE_SWITCH || defined LJ_POT_SWITCH
 576     /* We switch the LJ force */
 577     r_S0        = gmx_simd_mul_r(rsq_S0, rinv_S0);
 578     rsw_S0      = gmx_simd_max_r(gmx_simd_sub_r(r_S0, rswitch_S), zero_S);
 579     rsw2_S0     = gmx_simd_mul_r(rsw_S0, rsw_S0);
 580     rsw2_r_S0   = gmx_simd_mul_r(rsw2_S0, r_S0);
 581 #ifndef HALF_LJ
 582     r_S2        = gmx_simd_mul_r(rsq_S2, rinv_S2);
 583     rsw_S2      = gmx_simd_max_r(gmx_simd_sub_r(r_S2, rswitch_S), zero_S);
 584     rsw2_S2     = gmx_simd_mul_r(rsw_S2, rsw_S2);
 585     rsw2_r_S2   = gmx_simd_mul_r(rsw2_S2, r_S2);
 586 #endif
 587 #endif
 588
 589 #ifdef LJ_FORCE_SWITCH
 590
 591 #define add_fr_switch(fr, rsw, rsw2_r, c2, c3) gmx_simd_fmadd_r(gmx_simd_fmadd_r(c3, rsw, c2), rsw2_r, fr)
 592
 593     FrLJ6_S0    = gmx_simd_mul_r(c6_S0, add_fr_switch(rinvsix_S0, rsw_S0, rsw2_r_S0, p6_fc2_S, p6_fc3_S));
 594 #ifndef HALF_LJ
 595     FrLJ6_S2    = gmx_simd_mul_r(c6_S2, add_fr_switch(rinvsix_S2, rsw_S2, rsw2_r_S2, p6_fc2_S, p6_fc3_S));
 596 #endif
 597     FrLJ12_S0   = gmx_simd_mul_r(c12_S0, add_fr_switch(gmx_simd_mul_r(rinvsix_S0, rinvsix_S0), rsw_S0, rsw2_r_S0, p12_fc2_S, p12_fc3_S));
 598 #ifndef HALF_LJ
 599     FrLJ12_S2   = gmx_simd_mul_r(c12_S2, add_fr_switch(gmx_simd_mul_r(rinvsix_S2, rinvsix_S2), rsw_S2, rsw2_r_S2, p12_fc2_S, p12_fc3_S));
 600 #endif
 601 #undef add_fr_switch
 602 #endif /* LJ_FORCE_SWITCH */
 603
 604 #endif /* not LJ_COMB_LB */
 605
 606 #ifdef LJ_COMB_LB
 607     sir_S0      = gmx_simd_mul_r(sig_S0, rinv_S0);
 608 #ifndef HALF_LJ
 609     sir_S2      = gmx_simd_mul_r(sig_S2, rinv_S2);
 610 #endif
 611     sir2_S0     = gmx_simd_mul_r(sir_S0, sir_S0);
 612 #ifndef HALF_LJ
 613     sir2_S2     = gmx_simd_mul_r(sir_S2, sir_S2);
 614 #endif
 615     sir6_S0     = gmx_simd_mul_r(sir2_S0, gmx_simd_mul_r(sir2_S0, sir2_S0));
 616 #ifdef EXCL_FORCES
 617     sir6_S0     = gmx_simd_blendzero_r(sir6_S0, interact_S0);
 618 #endif
 619 #ifndef HALF_LJ
 620     sir6_S2     = gmx_simd_mul_r(sir2_S2, gmx_simd_mul_r(sir2_S2, sir2_S2));
 621 #ifdef EXCL_FORCES
 622     sir6_S2     = gmx_simd_blendzero_r(sir6_S2, interact_S2);
 623 #endif
 624 #endif
 625 #ifdef VDW_CUTOFF_CHECK
 626     sir6_S0     = gmx_simd_blendzero_r(sir6_S0, wco_vdw_S0);
 627 #ifndef HALF_LJ
 628     sir6_S2     = gmx_simd_blendzero_r(sir6_S2, wco_vdw_S2);
 629 #endif
 630 #endif
 631     FrLJ6_S0    = gmx_simd_mul_r(eps_S0, sir6_S0);
 632 #ifndef HALF_LJ
 633     FrLJ6_S2    = gmx_simd_mul_r(eps_S2, sir6_S2);
 634 #endif
 635     FrLJ12_S0   = gmx_simd_mul_r(FrLJ6_S0, sir6_S0);
 636 #ifndef HALF_LJ
 637     FrLJ12_S2   = gmx_simd_mul_r(FrLJ6_S2, sir6_S2);
 638 #endif
 639 #if defined CALC_ENERGIES
 640     /* We need C6 and C12 to calculate the LJ potential shift */
 641     sig2_S0     = gmx_simd_mul_r(sig_S0, sig_S0);
 642 #ifndef HALF_LJ
 643     sig2_S2     = gmx_simd_mul_r(sig_S2, sig_S2);
 644 #endif
 645     sig6_S0     = gmx_simd_mul_r(sig2_S0, gmx_simd_mul_r(sig2_S0, sig2_S0));
 646 #ifndef HALF_LJ
 647     sig6_S2     = gmx_simd_mul_r(sig2_S2, gmx_simd_mul_r(sig2_S2, sig2_S2));
 648 #endif
 649     c6_S0       = gmx_simd_mul_r(eps_S0, sig6_S0);
 650 #ifndef HALF_LJ
 651     c6_S2       = gmx_simd_mul_r(eps_S2, sig6_S2);
 652 #endif
 653     c12_S0      = gmx_simd_mul_r(c6_S0, sig6_S0);
 654 #ifndef HALF_LJ
 655     c12_S2      = gmx_simd_mul_r(c6_S2, sig6_S2);
 656 #endif
 657 #endif
 658 #endif /* LJ_COMB_LB */
 659
 660     /* Determine the total scalar LJ force*r */
 661     frLJ_S0     = gmx_simd_sub_r(FrLJ12_S0, FrLJ6_S0);
 662 #ifndef HALF_LJ
 663     frLJ_S2     = gmx_simd_sub_r(FrLJ12_S2, FrLJ6_S2);
 664 #endif
 665
 666 #if (defined LJ_CUT || defined LJ_FORCE_SWITCH) && defined CALC_ENERGIES
 667
 668 #ifdef LJ_CUT
 669     /* Calculate the LJ energies, with constant potential shift */
 670     VLJ6_S0     = gmx_simd_mul_r(sixth_S, gmx_simd_fmadd_r(c6_S0, p6_cpot_S, FrLJ6_S0));
 671 #ifndef HALF_LJ
 672     VLJ6_S2     = gmx_simd_mul_r(sixth_S, gmx_simd_fmadd_r(c6_S2, p6_cpot_S, FrLJ6_S2));
 673 #endif
 674     VLJ12_S0    = gmx_simd_mul_r(twelveth_S, gmx_simd_fmadd_r(c12_S0, p12_cpot_S, FrLJ12_S0));
 675 #ifndef HALF_LJ
 676     VLJ12_S2    = gmx_simd_mul_r(twelveth_S, gmx_simd_fmadd_r(c12_S2, p12_cpot_S, FrLJ12_S2));
 677 #endif
 678 #endif /* LJ_CUT */
 679
 680 #ifdef LJ_FORCE_SWITCH
 681 #define v_fswitch_pr(rsw, rsw2, c0, c3, c4) gmx_simd_fmadd_r(gmx_simd_fmadd_r(c4, rsw, c3), gmx_simd_mul_r(rsw2, rsw), c0)
 682
 683     VLJ6_S0     = gmx_simd_mul_r(c6_S0, gmx_simd_fmadd_r(sixth_S, rinvsix_S0, v_fswitch_pr(rsw_S0, rsw2_S0, p6_6cpot_S, p6_vc3_S, p6_vc4_S)));
 684 #ifndef HALF_LJ
 685     VLJ6_S2     = gmx_simd_mul_r(c6_S2, gmx_simd_fmadd_r(sixth_S, rinvsix_S2, v_fswitch_pr(rsw_S2, rsw2_S2, p6_6cpot_S, p6_vc3_S, p6_vc4_S)));
 686 #endif
 687     VLJ12_S0    = gmx_simd_mul_r(c12_S0, gmx_simd_fmadd_r(twelveth_S, gmx_simd_mul_r(rinvsix_S0, rinvsix_S0), v_fswitch_pr(rsw_S0, rsw2_S0, p12_12cpot_S, p12_vc3_S, p12_vc4_S)));
 688 #ifndef HALF_LJ
 689     VLJ12_S2    = gmx_simd_mul_r(c12_S2, gmx_simd_fmadd_r(twelveth_S, gmx_simd_mul_r(rinvsix_S2, rinvsix_S2), v_fswitch_pr(rsw_S2, rsw2_S2, p12_12cpot_S, p12_vc3_S, p12_vc4_S)));
 690 #endif
 691 #undef v_fswitch_pr
 692 #endif /* LJ_FORCE_SWITCH */
 693
 694     /* Add up the repulsion and dispersion */
 695     VLJ_S0      = gmx_simd_sub_r(VLJ12_S0, VLJ6_S0);
 696 #ifndef HALF_LJ
 697     VLJ_S2      = gmx_simd_sub_r(VLJ12_S2, VLJ6_S2);
 698 #endif
 699
 700 #endif /* (LJ_CUT || LJ_FORCE_SWITCH) && CALC_ENERGIES */
 701
 702 #ifdef LJ_POT_SWITCH
 703     /* We always need the potential, since it is needed for the force */
 704     VLJ_S0 = gmx_simd_fnmadd_r(sixth_S, FrLJ6_S0, gmx_simd_mul_r(twelveth_S, FrLJ12_S0));
 705 #ifndef HALF_LJ
 706     VLJ_S2 = gmx_simd_fnmadd_r(sixth_S, FrLJ6_S2, gmx_simd_mul_r(twelveth_S, FrLJ12_S2));
 707 #endif
 708
 709     {
 710         gmx_simd_real_t sw_S0, dsw_S0;
 711 #ifndef HALF_LJ
 712         gmx_simd_real_t sw_S2, dsw_S2;
 713 #endif
 714
 715 #define switch_pr(rsw, rsw2, c3, c4, c5) gmx_simd_fmadd_r(gmx_simd_fmadd_r(gmx_simd_fmadd_r(c5, rsw, c4), rsw, c3), gmx_simd_mul_r(rsw2, rsw), one_S)
 716 #define dswitch_pr(rsw, rsw2, c2, c3, c4) gmx_simd_mul_r(gmx_simd_fmadd_r(gmx_simd_fmadd_r(c4, rsw, c3), rsw, c2), rsw2)
 717
 718         sw_S0  = switch_pr(rsw_S0, rsw2_S0, swV3_S, swV4_S, swV5_S);
 719         dsw_S0 = dswitch_pr(rsw_S0, rsw2_S0, swF2_S, swF3_S, swF4_S);
 720 #ifndef HALF_LJ
 721         sw_S2  = switch_pr(rsw_S2, rsw2_S2, swV3_S, swV4_S, swV5_S);
 722         dsw_S2 = dswitch_pr(rsw_S2, rsw2_S2, swF2_S, swF3_S, swF4_S);
 723 #endif
 724         frLJ_S0 = gmx_simd_fnmadd_r(gmx_simd_mul_r(dsw_S0, VLJ_S0), r_S0, gmx_simd_mul_r(sw_S0, frLJ_S0));
 725 #ifndef HALF_LJ
 726         frLJ_S2 = gmx_simd_fnmadd_r(gmx_simd_mul_r(dsw_S2, VLJ_S2), r_S2, gmx_simd_mul_r(sw_S2, frLJ_S2));
 727 #endif
 728 #ifdef CALC_ENERGIES
 729         VLJ_S0  = gmx_simd_mul_r(sw_S0, VLJ_S0);
 730 #ifndef HALF_LJ
 731         VLJ_S2  = gmx_simd_mul_r(sw_S2, VLJ_S2);
 732 #endif
 733 #endif
 734
 735 #undef switch_pr
 736 #undef dswitch_pr
 737     }
 738 #endif /* LJ_POT_SWITCH */
 739
 740 #if defined CALC_ENERGIES && defined CHECK_EXCLS
 741     /* The potential shift should be removed for excluded pairs */
 742     VLJ_S0      = gmx_simd_blendzero_r(VLJ_S0, interact_S0);
 743 #ifndef HALF_LJ
 744     VLJ_S2      = gmx_simd_blendzero_r(VLJ_S2, interact_S2);
 745 #endif
 746 #endif
 747
 748 #ifdef LJ_EWALD_GEOM
 749     {
 750         gmx_simd_real_t c6s_j_S;
 751         gmx_simd_real_t c6grid_S0, rinvsix_nm_S0, cr2_S0, expmcr2_S0, poly_S0;
 752 #ifndef HALF_LJ
 753         gmx_simd_real_t c6grid_S2, rinvsix_nm_S2, cr2_S2, expmcr2_S2, poly_S2;
 754 #endif
 755 #ifdef CALC_ENERGIES
 756         gmx_simd_real_t sh_mask_S0;
 757 #ifndef HALF_LJ
 758         gmx_simd_real_t sh_mask_S2;
 759 #endif
 760 #endif
 761
 762         /* Determine C6 for the grid using the geometric combination rule */
 763         gmx_loaddh_pr(&c6s_j_S,  ljc+aj2+0);
 764         c6grid_S0       = gmx_simd_mul_r(c6s_S0, c6s_j_S);
 765 #ifndef HALF_LJ
 766         c6grid_S2       = gmx_simd_mul_r(c6s_S2, c6s_j_S);
 767 #endif
 768
 769 #ifdef CHECK_EXCLS
 770         /* Recalculate rinvsix without exclusion mask (compiler might optimize) */
 771         rinvsix_nm_S0 = gmx_simd_mul_r(rinvsq_S0, gmx_simd_mul_r(rinvsq_S0, rinvsq_S0));
 772 #ifndef HALF_LJ
 773         rinvsix_nm_S2 = gmx_simd_mul_r(rinvsq_S2, gmx_simd_mul_r(rinvsq_S2, rinvsq_S2));
 774 #endif
 775 #else
 776         /* We didn't use a mask, so we can copy */
 777         rinvsix_nm_S0 = rinvsix_S0;
 778 #ifndef HALF_LJ
 779         rinvsix_nm_S2 = rinvsix_S2;
 780 #endif
 781 #endif
 782
 783         cr2_S0        = gmx_simd_mul_r(lje_c2_S, rsq_S0);
 784 #ifndef HALF_LJ
 785         cr2_S2        = gmx_simd_mul_r(lje_c2_S, rsq_S2);
 786 #endif
 787         expmcr2_S0    = gmx_simd_exp_r(gmx_simd_mul_r(mone_S, cr2_S0));
 788 #ifndef HALF_LJ
 789         expmcr2_S2    = gmx_simd_exp_r(gmx_simd_mul_r(mone_S, cr2_S2));
 790 #endif
 791
 792         /* 1 + cr2 + 1/2*cr2^2 */
 793         poly_S0       = gmx_simd_fmadd_r(gmx_simd_fmadd_r(half_S, cr2_S0, one_S), cr2_S0, one_S);
 794 #ifndef HALF_LJ
 795         poly_S2       = gmx_simd_fmadd_r(gmx_simd_fmadd_r(half_S, cr2_S2, one_S), cr2_S2, one_S);
 796 #endif
 797
 798         /* We calculate LJ F*r = (6*C6)*(r^-6 - F_mesh/6), we use:
 799          * r^-6*cexp*(1 + cr2 + cr2^2/2 + cr2^3/6) = cexp*(r^-6*poly + c^6/6)
 800          */
 801         frLJ_S0       = gmx_simd_fmadd_r(c6grid_S0, gmx_simd_fnmadd_r(expmcr2_S0, gmx_simd_fmadd_r(rinvsix_nm_S0, poly_S0, lje_c6_6_S), rinvsix_nm_S0), frLJ_S0);
 802 #ifndef HALF_LJ
 803         frLJ_S2       = gmx_simd_fmadd_r(c6grid_S2, gmx_simd_fnmadd_r(expmcr2_S2, gmx_simd_fmadd_r(rinvsix_nm_S2, poly_S2, lje_c6_6_S), rinvsix_nm_S2), frLJ_S2);
 804 #endif
 805
 806 #ifdef CALC_ENERGIES
 807 #ifdef CHECK_EXCLS
 808         sh_mask_S0    = gmx_simd_blendzero_r(lje_vc_S, interact_S0);
 809 #ifndef HALF_LJ
 810         sh_mask_S2    = gmx_simd_blendzero_r(lje_vc_S, interact_S2);
 811 #endif
 812 #else
 813         sh_mask_S0    = lje_vc_S;
 814 #ifndef HALF_LJ
 815         sh_mask_S2    = lje_vc_S;
 816 #endif
 817 #endif
 818
 819         VLJ_S0        = gmx_simd_fmadd_r(gmx_simd_mul_r(sixth_S, c6grid_S0), gmx_simd_fmadd_r(rinvsix_nm_S0, gmx_simd_fnmadd_r(expmcr2_S0, poly_S0, one_S), sh_mask_S0), VLJ_S0);
 820 #ifndef HALF_LJ
 821         VLJ_S2        = gmx_simd_fmadd_r(gmx_simd_mul_r(sixth_S, c6grid_S2), gmx_simd_fmadd_r(rinvsix_nm_S2, gmx_simd_fnmadd_r(expmcr2_S2, poly_S2, one_S), sh_mask_S2), VLJ_S2);
 822 #endif
 823 #endif /* CALC_ENERGIES */
 824     }
 825 #endif /* LJ_EWALD_GEOM */
 826
 827 #if defined VDW_CUTOFF_CHECK
 828     /* frLJ is multiplied later by rinvsq, which is masked for the Coulomb
 829      * cut-off, but if the VdW cut-off is shorter, we need to mask with that.
 830      */
 831     frLJ_S0     = gmx_simd_blendzero_r(frLJ_S0, wco_vdw_S0);
 832 #ifndef HALF_LJ
 833     frLJ_S2     = gmx_simd_blendzero_r(frLJ_S2, wco_vdw_S2);
 834 #endif
 835 #endif
 836
 837 #ifdef CALC_ENERGIES
 838     /* The potential shift should be removed for pairs beyond cut-off */
 839     VLJ_S0      = gmx_simd_blendzero_r(VLJ_S0, wco_vdw_S0);
 840 #ifndef HALF_LJ
 841     VLJ_S2      = gmx_simd_blendzero_r(VLJ_S2, wco_vdw_S2);
 842 #endif
 843 #endif
 844
 845 #endif /* CALC_LJ */
 846
 847 #ifdef CALC_ENERGIES
 848 #ifdef ENERGY_GROUPS
 849     /* Extract the group pair index per j pair.
 850      * Energy groups are stored per i-cluster, so things get
 851      * complicated when the i- and j-cluster size don't match.
 852      */
 853     {
 854         int egps_j;
 855 #if UNROLLJ == 2
 856         egps_j    = nbat->energrp[cj>>1];
 857         egp_jj[0] = ((egps_j >> ((cj & 1)*egps_jshift)) & egps_jmask)*egps_jstride;
 858 #else
 859         /* We assume UNROLLI <= UNROLLJ */
 860         int jdi;
 861         for (jdi = 0; jdi < UNROLLJ/UNROLLI; jdi++)
 862         {
 863             int jj;
 864             egps_j = nbat->energrp[cj*(UNROLLJ/UNROLLI)+jdi];
 865             for (jj = 0; jj < (UNROLLI/2); jj++)
 866             {
 867                 egp_jj[jdi*(UNROLLI/2)+jj] = ((egps_j >> (jj*egps_jshift)) & egps_jmask)*egps_jstride;
 868             }
 869         }
 870 #endif
 871     }
 872 #endif
 873
 874 #ifdef CALC_COULOMB
 875 #ifndef ENERGY_GROUPS
 876     vctot_S      = gmx_simd_add_r(vctot_S, gmx_simd_add_r(vcoul_S0, vcoul_S2));
 877 #else
 878     add_ener_grp_halves(vcoul_S0, vctp[0], vctp[1], egp_jj);
 879     add_ener_grp_halves(vcoul_S2, vctp[2], vctp[3], egp_jj);
 880 #endif
 881 #endif
 882
 883 #ifdef CALC_LJ
 884 #ifndef ENERGY_GROUPS
 885     Vvdwtot_S    = gmx_simd_add_r(Vvdwtot_S,
 886 #ifndef HALF_LJ
 887                                   gmx_simd_add_r(VLJ_S0, VLJ_S2)
 888 #else
 889                                   VLJ_S0
 890 #endif
 891                                   );
 892 #else
 893     add_ener_grp_halves(VLJ_S0, vvdwtp[0], vvdwtp[1], egp_jj);
 894 #ifndef HALF_LJ
 895     add_ener_grp_halves(VLJ_S2, vvdwtp[2], vvdwtp[3], egp_jj);
 896 #endif
 897 #endif
 898 #endif /* CALC_LJ */
 899 #endif /* CALC_ENERGIES */
 900
 901 #ifdef CALC_LJ
 902 #ifdef CALC_COULOMB
 903     fscal_S0    = gmx_simd_mul_r(rinvsq_S0, gmx_simd_add_r(frcoul_S0, frLJ_S0));
 904 #else
 905     fscal_S0    = gmx_simd_mul_r(rinvsq_S0, frLJ_S0);
 906 #endif
 907 #else
 908     fscal_S0    = gmx_simd_mul_r(rinvsq_S0, frcoul_S0);
 909 #endif /* CALC_LJ */
 910 #if defined CALC_LJ && !defined HALF_LJ
 911 #ifdef CALC_COULOMB
 912     fscal_S2    = gmx_simd_mul_r(rinvsq_S2, gmx_simd_add_r(frcoul_S2, frLJ_S2));
 913 #else
 914     fscal_S2    = gmx_simd_mul_r(rinvsq_S2, frLJ_S2);
 915 #endif
 916 #else
 917     /* Atom 2 and 3 don't have LJ, so only add Coulomb forces */
 918     fscal_S2    = gmx_simd_mul_r(rinvsq_S2, frcoul_S2);
 919 #endif
 920
 921     /* Calculate temporary vectorial force */
 922     tx_S0       = gmx_simd_mul_r(fscal_S0, dx_S0);
 923     tx_S2       = gmx_simd_mul_r(fscal_S2, dx_S2);
 924     ty_S0       = gmx_simd_mul_r(fscal_S0, dy_S0);
 925     ty_S2       = gmx_simd_mul_r(fscal_S2, dy_S2);
 926     tz_S0       = gmx_simd_mul_r(fscal_S0, dz_S0);
 927     tz_S2       = gmx_simd_mul_r(fscal_S2, dz_S2);
 928
 929     /* Increment i atom force */
 930     fix_S0      = gmx_simd_add_r(fix_S0, tx_S0);
 931     fix_S2      = gmx_simd_add_r(fix_S2, tx_S2);
 932     fiy_S0      = gmx_simd_add_r(fiy_S0, ty_S0);
 933     fiy_S2      = gmx_simd_add_r(fiy_S2, ty_S2);
 934     fiz_S0      = gmx_simd_add_r(fiz_S0, tz_S0);
 935     fiz_S2      = gmx_simd_add_r(fiz_S2, tz_S2);
 936
 937     /* Decrement j atom force */
 938     gmx_load_hpr(&fjx_S, f+ajx);
 939     gmx_load_hpr(&fjy_S, f+ajy);
 940     gmx_load_hpr(&fjz_S, f+ajz);
 941     gmx_store_hpr(f+ajx, gmx_sub_hpr(fjx_S, gmx_sum4_hpr(tx_S0, tx_S2)));
 942     gmx_store_hpr(f+ajy, gmx_sub_hpr(fjy_S, gmx_sum4_hpr(ty_S0, ty_S2)));
 943     gmx_store_hpr(f+ajz, gmx_sub_hpr(fjz_S, gmx_sum4_hpr(tz_S0, tz_S2)));
 944 }
 945
 946 #undef  rinv_ex_S0
 947 #undef  rinv_ex_S2
 948
 949 #undef  wco_vdw_S0
 950 #undef  wco_vdw_S2
 951
 952 #undef  NBNXN_CUTOFF_USE_BLENDV
 953
 954 #undef  EXCL_FORCES