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