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
  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
  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 CUTOFF_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 CUTOFF_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 CUTOFF_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 #ifdef CALC_COULOMB
 118 #ifdef CHECK_EXCLS
 119     /* 1/r masked with the interaction mask */
 120     gmx_simd_real_t  rinv_ex_S0;
 121     gmx_simd_real_t  rinv_ex_S2;
 122 #endif
 123     gmx_simd_real_t  jq_S;
 124     gmx_simd_real_t  qq_S0;
 125     gmx_simd_real_t  qq_S2;
 126 #ifdef CALC_COUL_TAB
 127     /* The force (PME mesh force) we need to subtract from 1/r^2 */
 128     gmx_simd_real_t  fsub_S0;
 129     gmx_simd_real_t  fsub_S2;
 130 #endif
 131 #ifdef CALC_COUL_EWALD
 132     gmx_simd_real_t  brsq_S0, brsq_S2;
 133     gmx_simd_real_t  ewcorr_S0, ewcorr_S2;
 134 #endif
 135
 136     /* frcoul = (1/r - fsub)*r */
 137     gmx_simd_real_t  frcoul_S0;
 138     gmx_simd_real_t  frcoul_S2;
 139 #ifdef CALC_COUL_TAB
 140     /* For tables: r, rs=r/sp, rf=floor(rs), frac=rs-rf */
 141     gmx_simd_real_t         r_S0, rs_S0, rf_S0, frac_S0;
 142     gmx_simd_real_t         r_S2, rs_S2, rf_S2, frac_S2;
 143     /* Table index: rs truncated to an int */
 144     gmx_simd_int32_t        ti_S0, ti_S2;
 145     /* Linear force table values */
 146     gmx_simd_real_t         ctab0_S0, ctab1_S0;
 147     gmx_simd_real_t         ctab0_S2, ctab1_S2;
 148 #ifdef CALC_ENERGIES
 149     /* Quadratic energy table value */
 150     gmx_simd_real_t  ctabv_S0;
 151     gmx_simd_real_t  ctabv_S2;
 152 #endif
 153 #endif
 154 #if defined CALC_ENERGIES && (defined CALC_COUL_EWALD || defined CALC_COUL_TAB)
 155     /* The potential (PME mesh) we need to subtract from 1/r */
 156     gmx_simd_real_t  vc_sub_S0;
 157     gmx_simd_real_t  vc_sub_S2;
 158 #endif
 159 #ifdef CALC_ENERGIES
 160     /* Electrostatic potential */
 161     gmx_simd_real_t  vcoul_S0;
 162     gmx_simd_real_t  vcoul_S2;
 163 #endif
 164 #endif
 165     /* The force times 1/r */
 166     gmx_simd_real_t  fscal_S0;
 167     gmx_simd_real_t  fscal_S2;
 168
 169 #ifdef CALC_LJ
 170 #ifdef LJ_COMB_LB
 171     /* LJ sigma_j/2 and sqrt(epsilon_j) */
 172     gmx_simd_real_t  hsig_j_S, seps_j_S;
 173     /* LJ sigma_ij and epsilon_ij */
 174     gmx_simd_real_t  sig_S0, eps_S0;
 175 #ifndef HALF_LJ
 176     gmx_simd_real_t  sig_S2, eps_S2;
 177 #endif
 178 #ifdef CALC_ENERGIES
 179     gmx_simd_real_t  sig2_S0, sig6_S0;
 180 #ifndef HALF_LJ
 181     gmx_simd_real_t  sig2_S2, sig6_S2;
 182 #endif
 183 #endif /* LJ_COMB_LB */
 184 #endif /* CALC_LJ */
 185
 186 #ifdef LJ_COMB_GEOM
 187     gmx_simd_real_t  c6s_j_S, c12s_j_S;
 188 #endif
 189
 190 #if defined LJ_COMB_GEOM || defined LJ_COMB_LB
 191     /* Index for loading LJ parameters, complicated when interleaving */
 192     int         aj2;
 193 #endif
 194
 195 #ifndef FIX_LJ_C
 196     /* LJ C6 and C12 parameters, used with geometric comb. rule */
 197     gmx_simd_real_t  c6_S0, c12_S0;
 198 #ifndef HALF_LJ
 199     gmx_simd_real_t  c6_S2, c12_S2;
 200 #endif
 201 #endif
 202
 203     /* Intermediate variables for LJ calculation */
 204 #ifndef LJ_COMB_LB
 205     gmx_simd_real_t  rinvsix_S0;
 206 #ifndef HALF_LJ
 207     gmx_simd_real_t  rinvsix_S2;
 208 #endif
 209 #endif
 210 #ifdef LJ_COMB_LB
 211     gmx_simd_real_t  sir_S0, sir2_S0, sir6_S0;
 212 #ifndef HALF_LJ
 213     gmx_simd_real_t  sir_S2, sir2_S2, sir6_S2;
 214 #endif
 215 #endif
 216
 217     gmx_simd_real_t  FrLJ6_S0, FrLJ12_S0;
 218 #ifndef HALF_LJ
 219     gmx_simd_real_t  FrLJ6_S2, FrLJ12_S2;
 220 #endif
 221 #ifdef CALC_ENERGIES
 222     gmx_simd_real_t  VLJ6_S0, VLJ12_S0, VLJ_S0;
 223 #ifndef HALF_LJ
 224     gmx_simd_real_t  VLJ6_S2, VLJ12_S2, VLJ_S2;
 225 #endif
 226 #endif
 227 #endif /* CALC_LJ */
 228
 229     gmx_mm_hpr fjx_S, fjy_S, fjz_S;
 230
 231     /* j-cluster index */
 232     cj            = l_cj[cjind].cj;
 233
 234     /* Atom indices (of the first atom in the cluster) */
 235     aj            = cj*UNROLLJ;
 236 #if defined CALC_LJ && (defined LJ_COMB_GEOM || defined LJ_COMB_LB)
 237 #if UNROLLJ == STRIDE
 238     aj2           = aj*2;
 239 #else
 240     aj2           = (cj>>1)*2*STRIDE + (cj & 1)*UNROLLJ;
 241 #endif
 242 #endif
 243 #if UNROLLJ == STRIDE
 244     ajx           = aj*DIM;
 245 #else
 246     ajx           = (cj>>1)*DIM*STRIDE + (cj & 1)*UNROLLJ;
 247 #endif
 248     ajy           = ajx + STRIDE;
 249     ajz           = ajy + STRIDE;
 250
 251 #ifdef CHECK_EXCLS
 252     gmx_load_simd_2xnn_interactions(l_cj[cjind].excl,
 253                                     filter_S0, filter_S2,
 254                                     &interact_S0, &interact_S2);
 255 #endif /* CHECK_EXCLS */
 256
 257     /* load j atom coordinates */
 258     gmx_loaddh_pr(&jx_S, x+ajx);
 259     gmx_loaddh_pr(&jy_S, x+ajy);
 260     gmx_loaddh_pr(&jz_S, x+ajz);
 261
 262     /* Calculate distance */
 263     dx_S0       = gmx_simd_sub_r(ix_S0, jx_S);
 264     dy_S0       = gmx_simd_sub_r(iy_S0, jy_S);
 265     dz_S0       = gmx_simd_sub_r(iz_S0, jz_S);
 266     dx_S2       = gmx_simd_sub_r(ix_S2, jx_S);
 267     dy_S2       = gmx_simd_sub_r(iy_S2, jy_S);
 268     dz_S2       = gmx_simd_sub_r(iz_S2, jz_S);
 269
 270     /* rsq = dx*dx+dy*dy+dz*dz */
 271     rsq_S0      = gmx_simd_calc_rsq_r(dx_S0, dy_S0, dz_S0);
 272     rsq_S2      = gmx_simd_calc_rsq_r(dx_S2, dy_S2, dz_S2);
 273
 274 #ifndef CUTOFF_BLENDV
 275     wco_S0      = gmx_simd_cmplt_r(rsq_S0, rc2_S);
 276     wco_S2      = gmx_simd_cmplt_r(rsq_S2, rc2_S);
 277 #endif
 278
 279 #ifdef CHECK_EXCLS
 280 #ifdef EXCL_FORCES
 281     /* Only remove the (sub-)diagonal to avoid double counting */
 282 #if UNROLLJ == UNROLLI
 283     if (cj == ci_sh)
 284     {
 285         wco_S0  = gmx_simd_and_b(wco_S0, diagonal_mask_S0);
 286         wco_S2  = gmx_simd_and_b(wco_S2, diagonal_mask_S2);
 287     }
 288 #else
 289 #if UNROLLJ == 2*UNROLLI
 290     if (cj*2 == ci_sh)
 291     {
 292         wco_S0  = gmx_simd_and_b(wco_S0, diagonal_mask0_S0);
 293         wco_S2  = gmx_simd_and_b(wco_S2, diagonal_mask0_S2);
 294     }
 295     else if (cj*2 + 1 == ci_sh)
 296     {
 297         wco_S0  = gmx_simd_and_b(wco_S0, diagonal_mask1_S0);
 298         wco_S2  = gmx_simd_and_b(wco_S2, diagonal_mask1_S2);
 299     }
 300 #else
 301 #error "only UNROLLJ == UNROLLI*(1 or 2) currently supported in 2xnn kernels"
 302 #endif
 303 #endif
 304 #else /* EXCL_FORCES */
 305       /* No exclusion forces: remove all excluded atom pairs from the list */
 306     wco_S0      = gmx_simd_and_b(wco_S0, interact_S0);
 307     wco_S2      = gmx_simd_and_b(wco_S2, interact_S2);
 308 #endif
 309 #endif
 310
 311 #ifdef COUNT_PAIRS
 312     {
 313         int  i, j;
 314         real tmpa[2*GMX_SIMD_REAL_WIDTH], *tmp;
 315         tmp = gmx_simd_align_r(tmpa);
 316         for (i = 0; i < UNROLLI; i += 2)
 317         {
 318             gmx_simd_store_r(tmp, i == 0 ? wco_S0 : wco_S2);
 319             for (j = 0; j < 2*UNROLLJ; j++)
 320             {
 321                 if (!(tmp[j] == 0))
 322                 {
 323                     npair++;
 324                 }
 325             }
 326         }
 327     }
 328 #endif
 329
 330 #ifdef CHECK_EXCLS
 331     /* For excluded pairs add a small number to avoid r^-6 = NaN */
 332     rsq_S0      = gmx_masknot_add_pr(interact_S0, rsq_S0, avoid_sing_S);
 333     rsq_S2      = gmx_masknot_add_pr(interact_S2, rsq_S2, avoid_sing_S);
 334 #endif
 335
 336     /* Calculate 1/r */
 337     rinv_S0     = gmx_simd_invsqrt_r(rsq_S0);
 338     rinv_S2     = gmx_simd_invsqrt_r(rsq_S2);
 339
 340 #ifdef CALC_COULOMB
 341     /* Load parameters for j atom */
 342     gmx_loaddh_pr(&jq_S, q+aj);
 343     qq_S0       = gmx_simd_mul_r(iq_S0, jq_S);
 344     qq_S2       = gmx_simd_mul_r(iq_S2, jq_S);
 345 #endif
 346
 347 #ifdef CALC_LJ
 348
 349 #if !defined LJ_COMB_GEOM && !defined LJ_COMB_LB && !defined FIX_LJ_C
 350     load_lj_pair_params2(nbfp0, nbfp1, type, aj, &c6_S0, &c12_S0);
 351 #ifndef HALF_LJ
 352     load_lj_pair_params2(nbfp2, nbfp3, type, aj, &c6_S2, &c12_S2);
 353 #endif
 354 #endif /* not defined any LJ rule */
 355
 356 #ifdef LJ_COMB_GEOM
 357     gmx_loaddh_pr(&c6s_j_S,  ljc+aj2+0);
 358     gmx_loaddh_pr(&c12s_j_S, ljc+aj2+STRIDE);
 359     c6_S0       = gmx_simd_mul_r(c6s_S0, c6s_j_S );
 360 #ifndef HALF_LJ
 361     c6_S2       = gmx_simd_mul_r(c6s_S2, c6s_j_S );
 362 #endif
 363     c12_S0      = gmx_simd_mul_r(c12s_S0, c12s_j_S);
 364 #ifndef HALF_LJ
 365     c12_S2      = gmx_simd_mul_r(c12s_S2, c12s_j_S);
 366 #endif
 367 #endif /* LJ_COMB_GEOM */
 368
 369 #ifdef LJ_COMB_LB
 370     gmx_loaddh_pr(&hsig_j_S, ljc+aj2+0);
 371     gmx_loaddh_pr(&seps_j_S, ljc+aj2+STRIDE);
 372
 373     sig_S0      = gmx_simd_add_r(hsig_i_S0, hsig_j_S);
 374     eps_S0      = gmx_simd_mul_r(seps_i_S0, seps_j_S);
 375 #ifndef HALF_LJ
 376     sig_S2      = gmx_simd_add_r(hsig_i_S2, hsig_j_S);
 377     eps_S2      = gmx_simd_mul_r(seps_i_S2, seps_j_S);
 378 #endif
 379 #endif /* LJ_COMB_LB */
 380
 381 #endif /* CALC_LJ */
 382
 383 #ifndef CUTOFF_BLENDV
 384     rinv_S0     = gmx_simd_blendzero_r(rinv_S0, wco_S0);
 385     rinv_S2     = gmx_simd_blendzero_r(rinv_S2, wco_S2);
 386 #else
 387     /* We only need to mask for the cut-off: blendv is faster */
 388     rinv_S0     = gmx_simd_blendv_r(rinv_S0, zero_S, gmx_simd_sub_r(rc2_S, rsq_S0));
 389     rinv_S2     = gmx_simd_blendv_r(rinv_S2, zero_S, gmx_simd_sub_r(rc2_S, rsq_S2));
 390 #endif
 391
 392     rinvsq_S0   = gmx_simd_mul_r(rinv_S0, rinv_S0);
 393     rinvsq_S2   = gmx_simd_mul_r(rinv_S2, rinv_S2);
 394
 395 #ifdef CALC_COULOMB
 396     /* Note that here we calculate force*r, not the usual force/r.
 397      * This allows avoiding masking the reaction-field contribution,
 398      * as frcoul is later multiplied by rinvsq which has been
 399      * masked with the cut-off check.
 400      */
 401
 402 #ifdef EXCL_FORCES
 403     /* Only add 1/r for non-excluded atom pairs */
 404     rinv_ex_S0  = gmx_simd_blendzero_r(rinv_S0, interact_S0);
 405     rinv_ex_S2  = gmx_simd_blendzero_r(rinv_S2, interact_S2);
 406 #else
 407     /* No exclusion forces, we always need 1/r */
 408 #define     rinv_ex_S0    rinv_S0
 409 #define     rinv_ex_S2    rinv_S2
 410 #endif
 411
 412 #ifdef CALC_COUL_RF
 413     /* Electrostatic interactions */
 414     frcoul_S0   = gmx_simd_mul_r(qq_S0, gmx_simd_fmadd_r(rsq_S0, mrc_3_S, rinv_ex_S0));
 415     frcoul_S2   = gmx_simd_mul_r(qq_S2, gmx_simd_fmadd_r(rsq_S2, mrc_3_S, rinv_ex_S2));
 416
 417 #ifdef CALC_ENERGIES
 418     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)));
 419     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)));
 420 #endif
 421 #endif
 422
 423 #ifdef CALC_COUL_EWALD
 424     /* We need to mask (or limit) rsq for the cut-off,
 425      * as large distances can cause an overflow in gmx_pmecorrF/V.
 426      */
 427 #ifndef CUTOFF_BLENDV
 428     brsq_S0     = gmx_simd_mul_r(beta2_S, gmx_simd_blendzero_r(rsq_S0, wco_S0));
 429     brsq_S2     = gmx_simd_mul_r(beta2_S, gmx_simd_blendzero_r(rsq_S2, wco_S2));
 430 #else
 431     /* Strangely, putting mul on a separate line is slower (icc 13) */
 432     brsq_S0     = gmx_simd_mul_r(beta2_S, gmx_simd_blendv_r(rsq_S0, zero_S, gmx_simd_sub_r(rc2_S, rsq_S0)));
 433     brsq_S2     = gmx_simd_mul_r(beta2_S, gmx_simd_blendv_r(rsq_S2, zero_S, gmx_simd_sub_r(rc2_S, rsq_S2)));
 434 #endif
 435     ewcorr_S0   = gmx_simd_mul_r(gmx_simd_pmecorrF_r(brsq_S0), beta_S);
 436     ewcorr_S2   = gmx_simd_mul_r(gmx_simd_pmecorrF_r(brsq_S2), beta_S);
 437     frcoul_S0   = gmx_simd_mul_r(qq_S0, gmx_simd_fmadd_r(ewcorr_S0, brsq_S0, rinv_ex_S0));
 438     frcoul_S2   = gmx_simd_mul_r(qq_S2, gmx_simd_fmadd_r(ewcorr_S2, brsq_S2, rinv_ex_S2));
 439
 440 #ifdef CALC_ENERGIES
 441     vc_sub_S0   = gmx_simd_mul_r(gmx_simd_pmecorrV_r(brsq_S0), beta_S);
 442     vc_sub_S2   = gmx_simd_mul_r(gmx_simd_pmecorrV_r(brsq_S2), beta_S);
 443 #endif
 444
 445 #endif /* CALC_COUL_EWALD */
 446
 447 #ifdef CALC_COUL_TAB
 448     /* Electrostatic interactions */
 449     r_S0        = gmx_simd_mul_r(rsq_S0, rinv_S0);
 450     r_S2        = gmx_simd_mul_r(rsq_S2, rinv_S2);
 451     /* Convert r to scaled table units */
 452     rs_S0       = gmx_simd_mul_r(r_S0, invtsp_S);
 453     rs_S2       = gmx_simd_mul_r(r_S2, invtsp_S);
 454     /* Truncate scaled r to an int */
 455     ti_S0       = gmx_simd_cvtt_r2i(rs_S0);
 456     ti_S2       = gmx_simd_cvtt_r2i(rs_S2);
 457 #ifdef GMX_SIMD_HAVE_FLOOR
 458     rf_S0       = gmx_simd_floor_r(rs_S0);
 459     rf_S2       = gmx_simd_floor_r(rs_S2);
 460 #else
 461     rf_S0       = gmx_simd_cvt_i2r(ti_S0);
 462     rf_S2       = gmx_simd_cvt_i2r(ti_S2);
 463 #endif
 464     frac_S0     = gmx_simd_sub_r(rs_S0, rf_S0);
 465     frac_S2     = gmx_simd_sub_r(rs_S2, rf_S2);
 466
 467     /* Load and interpolate table forces and possibly energies.
 468      * Force and energy can be combined in one table, stride 4: FDV0
 469      * or in two separate tables with stride 1: F and V
 470      * Currently single precision uses FDV0, double F and V.
 471      */
 472 #ifndef CALC_ENERGIES
 473     load_table_f(tab_coul_F, ti_S0, ti0, &ctab0_S0, &ctab1_S0);
 474     load_table_f(tab_coul_F, ti_S2, ti2, &ctab0_S2, &ctab1_S2);
 475 #else
 476 #ifdef TAB_FDV0
 477     load_table_f_v(tab_coul_F, ti_S0, ti0, &ctab0_S0, &ctab1_S0, &ctabv_S0);
 478     load_table_f_v(tab_coul_F, ti_S2, ti2, &ctab0_S2, &ctab1_S2, &ctabv_S2);
 479 #else
 480     load_table_f_v(tab_coul_F, tab_coul_V, ti_S0, ti0, &ctab0_S0, &ctab1_S0, &ctabv_S0);
 481     load_table_f_v(tab_coul_F, tab_coul_V, ti_S2, ti2, &ctab0_S2, &ctab1_S2, &ctabv_S2);
 482 #endif
 483 #endif
 484     fsub_S0     = gmx_simd_add_r(ctab0_S0, gmx_simd_mul_r(frac_S0, ctab1_S0));
 485     fsub_S2     = gmx_simd_add_r(ctab0_S2, gmx_simd_mul_r(frac_S2, ctab1_S2));
 486     frcoul_S0   = gmx_simd_mul_r(qq_S0, gmx_simd_sub_r(rinv_ex_S0, gmx_simd_mul_r(fsub_S0, r_S0)));
 487     frcoul_S2   = gmx_simd_mul_r(qq_S2, gmx_simd_sub_r(rinv_ex_S2, gmx_simd_mul_r(fsub_S2, r_S2)));
 488
 489 #ifdef CALC_ENERGIES
 490     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)));
 491     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)));
 492 #endif
 493 #endif /* CALC_COUL_TAB */
 494
 495 #if defined CALC_ENERGIES && (defined CALC_COUL_EWALD || defined CALC_COUL_TAB)
 496 #ifndef NO_SHIFT_EWALD
 497     /* Add Ewald potential shift to vc_sub for convenience */
 498 #ifdef CHECK_EXCLS
 499     vc_sub_S0   = gmx_simd_add_r(vc_sub_S0, gmx_simd_blendzero_r(sh_ewald_S, interact_S0));
 500     vc_sub_S2   = gmx_simd_add_r(vc_sub_S2, gmx_simd_blendzero_r(sh_ewald_S, interact_S2));
 501 #else
 502     vc_sub_S0   = gmx_simd_add_r(vc_sub_S0, sh_ewald_S);
 503     vc_sub_S2   = gmx_simd_add_r(vc_sub_S2, sh_ewald_S);
 504 #endif
 505 #endif
 506
 507     vcoul_S0    = gmx_simd_mul_r(qq_S0, gmx_simd_sub_r(rinv_ex_S0, vc_sub_S0));
 508     vcoul_S2    = gmx_simd_mul_r(qq_S2, gmx_simd_sub_r(rinv_ex_S2, vc_sub_S2));
 509 #endif
 510
 511 #ifdef CALC_ENERGIES
 512     /* Mask energy for cut-off and diagonal */
 513     vcoul_S0    = gmx_simd_blendzero_r(vcoul_S0, wco_S0);
 514     vcoul_S2    = gmx_simd_blendzero_r(vcoul_S2, wco_S2);
 515 #endif
 516
 517 #endif /* CALC_COULOMB */
 518
 519 #ifdef CALC_LJ
 520     /* Lennard-Jones interaction */
 521
 522 #ifdef VDW_CUTOFF_CHECK
 523     wco_vdw_S0  = gmx_simd_cmplt_r(rsq_S0, rcvdw2_S);
 524 #ifndef HALF_LJ
 525     wco_vdw_S2  = gmx_simd_cmplt_r(rsq_S2, rcvdw2_S);
 526 #endif
 527 #else
 528     /* Same cut-off for Coulomb and VdW, reuse the registers */
 529 #define     wco_vdw_S0    wco_S0
 530 #define     wco_vdw_S2    wco_S2
 531 #endif
 532
 533 #ifndef LJ_COMB_LB
 534     rinvsix_S0  = gmx_simd_mul_r(rinvsq_S0, gmx_simd_mul_r(rinvsq_S0, rinvsq_S0));
 535 #ifdef EXCL_FORCES
 536     rinvsix_S0  = gmx_simd_blendzero_r(rinvsix_S0, interact_S0);
 537 #endif
 538 #ifndef HALF_LJ
 539     rinvsix_S2  = gmx_simd_mul_r(rinvsq_S2, gmx_simd_mul_r(rinvsq_S2, rinvsq_S2));
 540 #ifdef EXCL_FORCES
 541     rinvsix_S2  = gmx_simd_blendzero_r(rinvsix_S2, interact_S2);
 542 #endif
 543 #endif
 544 #ifdef VDW_CUTOFF_CHECK
 545     rinvsix_S0  = gmx_simd_blendzero_r(rinvsix_S0, wco_vdw_S0);
 546 #ifndef HALF_LJ
 547     rinvsix_S2  = gmx_simd_blendzero_r(rinvsix_S2, wco_vdw_S2);
 548 #endif
 549 #endif
 550     FrLJ6_S0    = gmx_simd_mul_r(c6_S0, rinvsix_S0);
 551 #ifndef HALF_LJ
 552     FrLJ6_S2    = gmx_simd_mul_r(c6_S2, rinvsix_S2);
 553 #endif
 554     FrLJ12_S0   = gmx_simd_mul_r(c12_S0, gmx_simd_mul_r(rinvsix_S0, rinvsix_S0));
 555 #ifndef HALF_LJ
 556     FrLJ12_S2   = gmx_simd_mul_r(c12_S2, gmx_simd_mul_r(rinvsix_S2, rinvsix_S2));
 557 #endif
 558 #endif /* not LJ_COMB_LB */
 559
 560 #ifdef LJ_COMB_LB
 561     sir_S0      = gmx_simd_mul_r(sig_S0, rinv_S0);
 562 #ifndef HALF_LJ
 563     sir_S2      = gmx_simd_mul_r(sig_S2, rinv_S2);
 564 #endif
 565     sir2_S0     = gmx_simd_mul_r(sir_S0, sir_S0);
 566 #ifndef HALF_LJ
 567     sir2_S2     = gmx_simd_mul_r(sir_S2, sir_S2);
 568 #endif
 569     sir6_S0     = gmx_simd_mul_r(sir2_S0, gmx_simd_mul_r(sir2_S0, sir2_S0));
 570 #ifdef EXCL_FORCES
 571     sir6_S0     = gmx_simd_blendzero_r(sir6_S0, interact_S0);
 572 #endif
 573 #ifndef HALF_LJ
 574     sir6_S2     = gmx_simd_mul_r(sir2_S2, gmx_simd_mul_r(sir2_S2, sir2_S2));
 575 #ifdef EXCL_FORCES
 576     sir6_S2     = gmx_simd_blendzero_r(sir6_S2, interact_S2);
 577 #endif
 578 #endif
 579 #ifdef VDW_CUTOFF_CHECK
 580     sir6_S0     = gmx_simd_blendzero_r(sir6_S0, wco_vdw_S0);
 581 #ifndef HALF_LJ
 582     sir6_S2     = gmx_simd_blendzero_r(sir6_S2, wco_vdw_S2);
 583 #endif
 584 #endif
 585     FrLJ6_S0    = gmx_simd_mul_r(eps_S0, sir6_S0);
 586 #ifndef HALF_LJ
 587     FrLJ6_S2    = gmx_simd_mul_r(eps_S2, sir6_S2);
 588 #endif
 589     FrLJ12_S0   = gmx_simd_mul_r(FrLJ6_S0, sir6_S0);
 590 #ifndef HALF_LJ
 591     FrLJ12_S2   = gmx_simd_mul_r(FrLJ6_S2, sir6_S2);
 592 #endif
 593 #if defined CALC_ENERGIES
 594     /* We need C6 and C12 to calculate the LJ potential shift */
 595     sig2_S0     = gmx_simd_mul_r(sig_S0, sig_S0);
 596 #ifndef HALF_LJ
 597     sig2_S2     = gmx_simd_mul_r(sig_S2, sig_S2);
 598 #endif
 599     sig6_S0     = gmx_simd_mul_r(sig2_S0, gmx_simd_mul_r(sig2_S0, sig2_S0));
 600 #ifndef HALF_LJ
 601     sig6_S2     = gmx_simd_mul_r(sig2_S2, gmx_simd_mul_r(sig2_S2, sig2_S2));
 602 #endif
 603     c6_S0       = gmx_simd_mul_r(eps_S0, sig6_S0);
 604 #ifndef HALF_LJ
 605     c6_S2       = gmx_simd_mul_r(eps_S2, sig6_S2);
 606 #endif
 607     c12_S0      = gmx_simd_mul_r(c6_S0, sig6_S0);
 608 #ifndef HALF_LJ
 609     c12_S2      = gmx_simd_mul_r(c6_S2, sig6_S2);
 610 #endif
 611 #endif
 612 #endif /* LJ_COMB_LB */
 613
 614 #endif /* CALC_LJ */
 615
 616 #ifdef CALC_ENERGIES
 617 #ifdef ENERGY_GROUPS
 618     /* Extract the group pair index per j pair.
 619      * Energy groups are stored per i-cluster, so things get
 620      * complicated when the i- and j-cluster size don't match.
 621      */
 622     {
 623         int egps_j;
 624 #if UNROLLJ == 2
 625         egps_j    = nbat->energrp[cj>>1];
 626         egp_jj[0] = ((egps_j >> ((cj & 1)*egps_jshift)) & egps_jmask)*egps_jstride;
 627 #else
 628         /* We assume UNROLLI <= UNROLLJ */
 629         int jdi;
 630         for (jdi = 0; jdi < UNROLLJ/UNROLLI; jdi++)
 631         {
 632             int jj;
 633             egps_j = nbat->energrp[cj*(UNROLLJ/UNROLLI)+jdi];
 634             for (jj = 0; jj < (UNROLLI/2); jj++)
 635             {
 636                 egp_jj[jdi*(UNROLLI/2)+jj] = ((egps_j >> (jj*egps_jshift)) & egps_jmask)*egps_jstride;
 637             }
 638         }
 639 #endif
 640     }
 641 #endif
 642
 643 #ifdef CALC_COULOMB
 644 #ifndef ENERGY_GROUPS
 645     vctot_S      = gmx_simd_add_r(vctot_S, gmx_simd_add_r(vcoul_S0, vcoul_S2));
 646 #else
 647     add_ener_grp_halves(vcoul_S0, vctp[0], vctp[1], egp_jj);
 648     add_ener_grp_halves(vcoul_S2, vctp[2], vctp[3], egp_jj);
 649 #endif
 650 #endif
 651
 652 #ifdef CALC_LJ
 653     /* Calculate the LJ energies */
 654     VLJ6_S0     = gmx_simd_mul_r(sixth_S, gmx_simd_sub_r(FrLJ6_S0, gmx_simd_mul_r(c6_S0, sh_invrc6_S)));
 655 #ifndef HALF_LJ
 656     VLJ6_S2     = gmx_simd_mul_r(sixth_S, gmx_simd_sub_r(FrLJ6_S2, gmx_simd_mul_r(c6_S2, sh_invrc6_S)));
 657 #endif
 658     VLJ12_S0    = gmx_simd_mul_r(twelveth_S, gmx_simd_sub_r(FrLJ12_S0, gmx_simd_mul_r(c12_S0, sh_invrc12_S)));
 659 #ifndef HALF_LJ
 660     VLJ12_S2    = gmx_simd_mul_r(twelveth_S, gmx_simd_sub_r(FrLJ12_S2, gmx_simd_mul_r(c12_S2, sh_invrc12_S)));
 661 #endif
 662
 663     VLJ_S0      = gmx_simd_sub_r(VLJ12_S0, VLJ6_S0);
 664 #ifndef HALF_LJ
 665     VLJ_S2      = gmx_simd_sub_r(VLJ12_S2, VLJ6_S2);
 666 #endif
 667     /* The potential shift should be removed for pairs beyond cut-off */
 668     VLJ_S0      = gmx_simd_blendzero_r(VLJ_S0, wco_vdw_S0);
 669 #ifndef HALF_LJ
 670     VLJ_S2      = gmx_simd_blendzero_r(VLJ_S2, wco_vdw_S2);
 671 #endif
 672 #ifdef CHECK_EXCLS
 673     /* The potential shift should be removed for excluded pairs */
 674     VLJ_S0      = gmx_simd_blendzero_r(VLJ_S0, interact_S0);
 675 #ifndef HALF_LJ
 676     VLJ_S2      = gmx_simd_blendzero_r(VLJ_S2, interact_S2);
 677 #endif
 678 #endif
 679 #ifndef ENERGY_GROUPS
 680     Vvdwtot_S    = gmx_simd_add_r(Vvdwtot_S,
 681 #ifndef HALF_LJ
 682                                   gmx_simd_add_r(VLJ_S0, VLJ_S2)
 683 #else
 684                                   VLJ_S0
 685 #endif
 686                                   );
 687 #else
 688     add_ener_grp_halves(VLJ_S0, vvdwtp[0], vvdwtp[1], egp_jj);
 689 #ifndef HALF_LJ
 690     add_ener_grp_halves(VLJ_S2, vvdwtp[2], vvdwtp[3], egp_jj);
 691 #endif
 692 #endif
 693 #endif /* CALC_LJ */
 694 #endif /* CALC_ENERGIES */
 695
 696 #ifdef CALC_LJ
 697 #ifdef CALC_COULOMB
 698     fscal_S0    = gmx_simd_mul_r(rinvsq_S0,
 699                                  gmx_simd_add_r(frcoul_S0,
 700                                                 gmx_simd_sub_r(FrLJ12_S0, FrLJ6_S0)));
 701 #else
 702     fscal_S0    = gmx_simd_mul_r(rinvsq_S0,
 703                                  (
 704                                      gmx_simd_sub_r(FrLJ12_S0, FrLJ6_S0)));
 705 #endif
 706 #else
 707     fscal_S0    = gmx_simd_mul_r(rinvsq_S0, frcoul_S0);
 708 #endif /* CALC_LJ */
 709 #if defined CALC_LJ && !defined HALF_LJ
 710 #ifdef CALC_COULOMB
 711     fscal_S2    = gmx_simd_mul_r(rinvsq_S2,
 712                                  gmx_simd_add_r(frcoul_S2,
 713                                                 gmx_simd_sub_r(FrLJ12_S2, FrLJ6_S2)));
 714 #else
 715     fscal_S2    = gmx_simd_mul_r(rinvsq_S2,
 716                                  (
 717                                      gmx_simd_sub_r(FrLJ12_S2, FrLJ6_S2)));
 718 #endif
 719 #else
 720     /* Atom 2 and 3 don't have LJ, so only add Coulomb forces */
 721     fscal_S2    = gmx_simd_mul_r(rinvsq_S2, frcoul_S2);
 722 #endif
 723
 724     /* Calculate temporary vectorial force */
 725     tx_S0       = gmx_simd_mul_r(fscal_S0, dx_S0);
 726     tx_S2       = gmx_simd_mul_r(fscal_S2, dx_S2);
 727     ty_S0       = gmx_simd_mul_r(fscal_S0, dy_S0);
 728     ty_S2       = gmx_simd_mul_r(fscal_S2, dy_S2);
 729     tz_S0       = gmx_simd_mul_r(fscal_S0, dz_S0);
 730     tz_S2       = gmx_simd_mul_r(fscal_S2, dz_S2);
 731
 732     /* Increment i atom force */
 733     fix_S0      = gmx_simd_add_r(fix_S0, tx_S0);
 734     fix_S2      = gmx_simd_add_r(fix_S2, tx_S2);
 735     fiy_S0      = gmx_simd_add_r(fiy_S0, ty_S0);
 736     fiy_S2      = gmx_simd_add_r(fiy_S2, ty_S2);
 737     fiz_S0      = gmx_simd_add_r(fiz_S0, tz_S0);
 738     fiz_S2      = gmx_simd_add_r(fiz_S2, tz_S2);
 739
 740     /* Decrement j atom force */
 741     gmx_load_hpr(&fjx_S, f+ajx);
 742     gmx_load_hpr(&fjy_S, f+ajy);
 743     gmx_load_hpr(&fjz_S, f+ajz);
 744     gmx_store_hpr(f+ajx, gmx_sub_hpr(fjx_S, gmx_sum4_hpr(tx_S0, tx_S2)));
 745     gmx_store_hpr(f+ajy, gmx_sub_hpr(fjy_S, gmx_sum4_hpr(ty_S0, ty_S2)));
 746     gmx_store_hpr(f+ajz, gmx_sub_hpr(fjz_S, gmx_sum4_hpr(tz_S0, tz_S2)));
 747 }
 748
 749 #undef  rinv_ex_S0
 750 #undef  rinv_ex_S2
 751
 752 #undef  wco_vdw_S0
 753 #undef  wco_vdw_S2
 754
 755 #undef  CUTOFF_BLENDV
 756
 757 #undef  EXCL_FORCES