src/gromacs/gmxlib/nonbonded/nb_kernel_avx_128_fma_double/nb_kernel_ElecEwSh_VdwLJEwSh_GeomW3P1_avx_128_fma_double.c

   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  * Note: this file was generated by the GROMACS avx_128_fma_double kernel generator.
  37  */
  38 #include "config.h"
  39
  40 #include <math.h>
  41
  42 #include "../nb_kernel.h"
  43 #include "gromacs/legacyheaders/types/simple.h"
  44 #include "gromacs/math/vec.h"
  45 #include "gromacs/legacyheaders/nrnb.h"
  46
  47 #include "gromacs/simd/math_x86_avx_128_fma_double.h"
  48 #include "kernelutil_x86_avx_128_fma_double.h"
  49
  50 /*
  51  * Gromacs nonbonded kernel:   nb_kernel_ElecEwSh_VdwLJEwSh_GeomW3P1_VF_avx_128_fma_double
  52  * Electrostatics interaction: Ewald
  53  * VdW interaction:            LJEwald
  54  * Geometry:                   Water3-Particle
  55  * Calculate force/pot:        PotentialAndForce
  56  */
  57 void
  58 nb_kernel_ElecEwSh_VdwLJEwSh_GeomW3P1_VF_avx_128_fma_double
  59                     (t_nblist                    * gmx_restrict       nlist,
  60                      rvec                        * gmx_restrict          xx,
  61                      rvec                        * gmx_restrict          ff,
  62                      t_forcerec                  * gmx_restrict          fr,
  63                      t_mdatoms                   * gmx_restrict     mdatoms,
  64                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
  65                      t_nrnb                      * gmx_restrict        nrnb)
  66 {
  67     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
  68      * just 0 for non-waters.
  69      * Suffixes A,B refer to j loop unrolling done with SSE double precision, e.g. for the two different
  70      * jnr indices corresponding to data put in the four positions in the SIMD register.
  71      */
  72     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
  73     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
  74     int              jnrA,jnrB;
  75     int              j_coord_offsetA,j_coord_offsetB;
  76     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
  77     real             rcutoff_scalar;
  78     real             *shiftvec,*fshift,*x,*f;
  79     __m128d          tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
  80     int              vdwioffset0;
  81     __m128d          ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
  82     int              vdwioffset1;
  83     __m128d          ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
  84     int              vdwioffset2;
  85     __m128d          ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
  86     int              vdwjidx0A,vdwjidx0B;
  87     __m128d          jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
  88     __m128d          dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
  89     __m128d          dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
  90     __m128d          dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
  91     __m128d          velec,felec,velecsum,facel,crf,krf,krf2;
  92     real             *charge;
  93     int              nvdwtype;
  94     __m128d          rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
  95     int              *vdwtype;
  96     real             *vdwparam;
  97     __m128d          one_sixth   = _mm_set1_pd(1.0/6.0);
  98     __m128d          one_twelfth = _mm_set1_pd(1.0/12.0);
  99     __m128d           c6grid_00;
 100     __m128d           c6grid_10;
 101     __m128d           c6grid_20;
 102     real             *vdwgridparam;
 103     __m128d           ewclj,ewclj2,ewclj6,ewcljrsq,poly,exponent,f6A,f6B,sh_lj_ewald;
 104     __m128d           one_half  = _mm_set1_pd(0.5);
 105     __m128d           minus_one = _mm_set1_pd(-1.0);
 106     __m128i          ewitab;
 107     __m128d          ewtabscale,eweps,twoeweps,sh_ewald,ewrt,ewtabhalfspace,ewtabF,ewtabFn,ewtabD,ewtabV;
 108     real             *ewtab;
 109     __m128d          dummy_mask,cutoff_mask;
 110     __m128d          signbit   = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
 111     __m128d          one     = _mm_set1_pd(1.0);
 112     __m128d          two     = _mm_set1_pd(2.0);
 113     x                = xx[0];
 114     f                = ff[0];
 115
 116     nri              = nlist->nri;
 117     iinr             = nlist->iinr;
 118     jindex           = nlist->jindex;
 119     jjnr             = nlist->jjnr;
 120     shiftidx         = nlist->shift;
 121     gid              = nlist->gid;
 122     shiftvec         = fr->shift_vec[0];
 123     fshift           = fr->fshift[0];
 124     facel            = _mm_set1_pd(fr->epsfac);
 125     charge           = mdatoms->chargeA;
 126     nvdwtype         = fr->ntype;
 127     vdwparam         = fr->nbfp;
 128     vdwtype          = mdatoms->typeA;
 129     vdwgridparam     = fr->ljpme_c6grid;
 130     sh_lj_ewald      = _mm_set1_pd(fr->ic->sh_lj_ewald);
 131     ewclj            = _mm_set1_pd(fr->ewaldcoeff_lj);
 132     ewclj2           = _mm_mul_pd(minus_one,_mm_mul_pd(ewclj,ewclj));
 133
 134     sh_ewald         = _mm_set1_pd(fr->ic->sh_ewald);
 135     ewtab            = fr->ic->tabq_coul_FDV0;
 136     ewtabscale       = _mm_set1_pd(fr->ic->tabq_scale);
 137     ewtabhalfspace   = _mm_set1_pd(0.5/fr->ic->tabq_scale);
 138
 139     /* Setup water-specific parameters */
 140     inr              = nlist->iinr[0];
 141     iq0              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+0]));
 142     iq1              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+1]));
 143     iq2              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+2]));
 144     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
 145
 146     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
 147     rcutoff_scalar   = fr->rcoulomb;
 148     rcutoff          = _mm_set1_pd(rcutoff_scalar);
 149     rcutoff2         = _mm_mul_pd(rcutoff,rcutoff);
 150
 151     sh_vdw_invrcut6  = _mm_set1_pd(fr->ic->sh_invrc6);
 152     rvdw             = _mm_set1_pd(fr->rvdw);
 153
 154     /* Avoid stupid compiler warnings */
 155     jnrA = jnrB = 0;
 156     j_coord_offsetA = 0;
 157     j_coord_offsetB = 0;
 158
 159     outeriter        = 0;
 160     inneriter        = 0;
 161
 162     /* Start outer loop over neighborlists */
 163     for(iidx=0; iidx<nri; iidx++)
 164     {
 165         /* Load shift vector for this list */
 166         i_shift_offset   = DIM*shiftidx[iidx];
 167
 168         /* Load limits for loop over neighbors */
 169         j_index_start    = jindex[iidx];
 170         j_index_end      = jindex[iidx+1];
 171
 172         /* Get outer coordinate index */
 173         inr              = iinr[iidx];
 174         i_coord_offset   = DIM*inr;
 175
 176         /* Load i particle coords and add shift vector */
 177         gmx_mm_load_shift_and_3rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,
 178                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2);
 179
 180         fix0             = _mm_setzero_pd();
 181         fiy0             = _mm_setzero_pd();
 182         fiz0             = _mm_setzero_pd();
 183         fix1             = _mm_setzero_pd();
 184         fiy1             = _mm_setzero_pd();
 185         fiz1             = _mm_setzero_pd();
 186         fix2             = _mm_setzero_pd();
 187         fiy2             = _mm_setzero_pd();
 188         fiz2             = _mm_setzero_pd();
 189
 190         /* Reset potential sums */
 191         velecsum         = _mm_setzero_pd();
 192         vvdwsum          = _mm_setzero_pd();
 193
 194         /* Start inner kernel loop */
 195         for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
 196         {
 197
 198             /* Get j neighbor index, and coordinate index */
 199             jnrA             = jjnr[jidx];
 200             jnrB             = jjnr[jidx+1];
 201             j_coord_offsetA  = DIM*jnrA;
 202             j_coord_offsetB  = DIM*jnrB;
 203
 204             /* load j atom coordinates */
 205             gmx_mm_load_1rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
 206                                               &jx0,&jy0,&jz0);
 207
 208             /* Calculate displacement vector */
 209             dx00             = _mm_sub_pd(ix0,jx0);
 210             dy00             = _mm_sub_pd(iy0,jy0);
 211             dz00             = _mm_sub_pd(iz0,jz0);
 212             dx10             = _mm_sub_pd(ix1,jx0);
 213             dy10             = _mm_sub_pd(iy1,jy0);
 214             dz10             = _mm_sub_pd(iz1,jz0);
 215             dx20             = _mm_sub_pd(ix2,jx0);
 216             dy20             = _mm_sub_pd(iy2,jy0);
 217             dz20             = _mm_sub_pd(iz2,jz0);
 218
 219             /* Calculate squared distance and things based on it */
 220             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
 221             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
 222             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
 223
 224             rinv00           = gmx_mm_invsqrt_pd(rsq00);
 225             rinv10           = gmx_mm_invsqrt_pd(rsq10);
 226             rinv20           = gmx_mm_invsqrt_pd(rsq20);
 227
 228             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
 229             rinvsq10         = _mm_mul_pd(rinv10,rinv10);
 230             rinvsq20         = _mm_mul_pd(rinv20,rinv20);
 231
 232             /* Load parameters for j particles */
 233             jq0              = gmx_mm_load_2real_swizzle_pd(charge+jnrA+0,charge+jnrB+0);
 234             vdwjidx0A        = 2*vdwtype[jnrA+0];
 235             vdwjidx0B        = 2*vdwtype[jnrB+0];
 236
 237             fjx0             = _mm_setzero_pd();
 238             fjy0             = _mm_setzero_pd();
 239             fjz0             = _mm_setzero_pd();
 240
 241             /**************************
 242              * CALCULATE INTERACTIONS *
 243              **************************/
 244
 245             if (gmx_mm_any_lt(rsq00,rcutoff2))
 246             {
 247
 248             r00              = _mm_mul_pd(rsq00,rinv00);
 249
 250             /* Compute parameters for interactions between i and j atoms */
 251             qq00             = _mm_mul_pd(iq0,jq0);
 252             gmx_mm_load_2pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,
 253                                          vdwparam+vdwioffset0+vdwjidx0B,&c6_00,&c12_00);
 254             c6grid_00       = gmx_mm_load_2real_swizzle_pd(vdwgridparam+vdwioffset0+vdwjidx0A,
 255                                                                vdwgridparam+vdwioffset0+vdwjidx0B);
 256
 257             /* EWALD ELECTROSTATICS */
 258
 259             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 260             ewrt             = _mm_mul_pd(r00,ewtabscale);
 261             ewitab           = _mm_cvttpd_epi32(ewrt);
 262 #ifdef __XOP__
 263             eweps            = _mm_frcz_pd(ewrt);
 264 #else
 265             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 266 #endif
 267             twoeweps         = _mm_add_pd(eweps,eweps);
 268             ewitab           = _mm_slli_epi32(ewitab,2);
 269             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 270             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 271             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 272             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 273             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
 274             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 275             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 276             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 277             velec            = _mm_mul_pd(qq00,_mm_sub_pd(_mm_sub_pd(rinv00,sh_ewald),velec));
 278             felec            = _mm_mul_pd(_mm_mul_pd(qq00,rinv00),_mm_sub_pd(rinvsq00,felec));
 279
 280             /* Analytical LJ-PME */
 281             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
 282             ewcljrsq         = _mm_mul_pd(ewclj2,rsq00);
 283             ewclj6           = _mm_mul_pd(ewclj2,_mm_mul_pd(ewclj2,ewclj2));
 284             exponent         = gmx_simd_exp_d(ewcljrsq);
 285             /* poly = exp(-(beta*r)^2) * (1 + (beta*r)^2 + (beta*r)^4 /2) */
 286             poly             = _mm_mul_pd(exponent,_mm_macc_pd(_mm_mul_pd(ewcljrsq,ewcljrsq),one_half,_mm_sub_pd(one,ewcljrsq)));
 287             /* vvdw6 = [C6 - C6grid * (1-poly)]/r6 */
 288             vvdw6            = _mm_mul_pd(_mm_macc_pd(-c6grid_00,_mm_sub_pd(one,poly),c6_00),rinvsix);
 289             vvdw12           = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
 290             vvdw             = _mm_msub_pd(_mm_nmacc_pd(c12_00,_mm_mul_pd(sh_vdw_invrcut6,sh_vdw_invrcut6),vvdw12),one_twelfth,
 291                                _mm_mul_pd(_mm_sub_pd(vvdw6,_mm_macc_pd(c6grid_00,sh_lj_ewald,_mm_mul_pd(c6_00,sh_vdw_invrcut6))),one_sixth));
 292             /* fvdw = vvdw12/r - (vvdw6/r + (C6grid * exponent * beta^6)/r) */
 293             fvdw             = _mm_mul_pd(_mm_add_pd(vvdw12,_mm_msub_pd(_mm_mul_pd(c6grid_00,one_sixth),_mm_mul_pd(exponent,ewclj6),vvdw6)),rinvsq00);
 294
 295             cutoff_mask      = _mm_cmplt_pd(rsq00,rcutoff2);
 296
 297             /* Update potential sum for this i atom from the interaction with this j atom. */
 298             velec            = _mm_and_pd(velec,cutoff_mask);
 299             velecsum         = _mm_add_pd(velecsum,velec);
 300             vvdw             = _mm_and_pd(vvdw,cutoff_mask);
 301             vvdwsum          = _mm_add_pd(vvdwsum,vvdw);
 302
 303             fscal            = _mm_add_pd(felec,fvdw);
 304
 305             fscal            = _mm_and_pd(fscal,cutoff_mask);
 306
 307             /* Update vectorial force */
 308             fix0             = _mm_macc_pd(dx00,fscal,fix0);
 309             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
 310             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
 311
 312             fjx0             = _mm_macc_pd(dx00,fscal,fjx0);
 313             fjy0             = _mm_macc_pd(dy00,fscal,fjy0);
 314             fjz0             = _mm_macc_pd(dz00,fscal,fjz0);
 315
 316             }
 317
 318             /**************************
 319              * CALCULATE INTERACTIONS *
 320              **************************/
 321
 322             if (gmx_mm_any_lt(rsq10,rcutoff2))
 323             {
 324
 325             r10              = _mm_mul_pd(rsq10,rinv10);
 326
 327             /* Compute parameters for interactions between i and j atoms */
 328             qq10             = _mm_mul_pd(iq1,jq0);
 329
 330             /* EWALD ELECTROSTATICS */
 331
 332             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 333             ewrt             = _mm_mul_pd(r10,ewtabscale);
 334             ewitab           = _mm_cvttpd_epi32(ewrt);
 335 #ifdef __XOP__
 336             eweps            = _mm_frcz_pd(ewrt);
 337 #else
 338             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 339 #endif
 340             twoeweps         = _mm_add_pd(eweps,eweps);
 341             ewitab           = _mm_slli_epi32(ewitab,2);
 342             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 343             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 344             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 345             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 346             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
 347             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 348             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 349             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 350             velec            = _mm_mul_pd(qq10,_mm_sub_pd(_mm_sub_pd(rinv10,sh_ewald),velec));
 351             felec            = _mm_mul_pd(_mm_mul_pd(qq10,rinv10),_mm_sub_pd(rinvsq10,felec));
 352
 353             cutoff_mask      = _mm_cmplt_pd(rsq10,rcutoff2);
 354
 355             /* Update potential sum for this i atom from the interaction with this j atom. */
 356             velec            = _mm_and_pd(velec,cutoff_mask);
 357             velecsum         = _mm_add_pd(velecsum,velec);
 358
 359             fscal            = felec;
 360
 361             fscal            = _mm_and_pd(fscal,cutoff_mask);
 362
 363             /* Update vectorial force */
 364             fix1             = _mm_macc_pd(dx10,fscal,fix1);
 365             fiy1             = _mm_macc_pd(dy10,fscal,fiy1);
 366             fiz1             = _mm_macc_pd(dz10,fscal,fiz1);
 367
 368             fjx0             = _mm_macc_pd(dx10,fscal,fjx0);
 369             fjy0             = _mm_macc_pd(dy10,fscal,fjy0);
 370             fjz0             = _mm_macc_pd(dz10,fscal,fjz0);
 371
 372             }
 373
 374             /**************************
 375              * CALCULATE INTERACTIONS *
 376              **************************/
 377
 378             if (gmx_mm_any_lt(rsq20,rcutoff2))
 379             {
 380
 381             r20              = _mm_mul_pd(rsq20,rinv20);
 382
 383             /* Compute parameters for interactions between i and j atoms */
 384             qq20             = _mm_mul_pd(iq2,jq0);
 385
 386             /* EWALD ELECTROSTATICS */
 387
 388             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 389             ewrt             = _mm_mul_pd(r20,ewtabscale);
 390             ewitab           = _mm_cvttpd_epi32(ewrt);
 391 #ifdef __XOP__
 392             eweps            = _mm_frcz_pd(ewrt);
 393 #else
 394             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 395 #endif
 396             twoeweps         = _mm_add_pd(eweps,eweps);
 397             ewitab           = _mm_slli_epi32(ewitab,2);
 398             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 399             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 400             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 401             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 402             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
 403             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 404             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 405             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 406             velec            = _mm_mul_pd(qq20,_mm_sub_pd(_mm_sub_pd(rinv20,sh_ewald),velec));
 407             felec            = _mm_mul_pd(_mm_mul_pd(qq20,rinv20),_mm_sub_pd(rinvsq20,felec));
 408
 409             cutoff_mask      = _mm_cmplt_pd(rsq20,rcutoff2);
 410
 411             /* Update potential sum for this i atom from the interaction with this j atom. */
 412             velec            = _mm_and_pd(velec,cutoff_mask);
 413             velecsum         = _mm_add_pd(velecsum,velec);
 414
 415             fscal            = felec;
 416
 417             fscal            = _mm_and_pd(fscal,cutoff_mask);
 418
 419             /* Update vectorial force */
 420             fix2             = _mm_macc_pd(dx20,fscal,fix2);
 421             fiy2             = _mm_macc_pd(dy20,fscal,fiy2);
 422             fiz2             = _mm_macc_pd(dz20,fscal,fiz2);
 423
 424             fjx0             = _mm_macc_pd(dx20,fscal,fjx0);
 425             fjy0             = _mm_macc_pd(dy20,fscal,fjy0);
 426             fjz0             = _mm_macc_pd(dz20,fscal,fjz0);
 427
 428             }
 429
 430             gmx_mm_decrement_1rvec_2ptr_swizzle_pd(f+j_coord_offsetA,f+j_coord_offsetB,fjx0,fjy0,fjz0);
 431
 432             /* Inner loop uses 179 flops */
 433         }
 434
 435         if(jidx<j_index_end)
 436         {
 437
 438             jnrA             = jjnr[jidx];
 439             j_coord_offsetA  = DIM*jnrA;
 440
 441             /* load j atom coordinates */
 442             gmx_mm_load_1rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
 443                                               &jx0,&jy0,&jz0);
 444
 445             /* Calculate displacement vector */
 446             dx00             = _mm_sub_pd(ix0,jx0);
 447             dy00             = _mm_sub_pd(iy0,jy0);
 448             dz00             = _mm_sub_pd(iz0,jz0);
 449             dx10             = _mm_sub_pd(ix1,jx0);
 450             dy10             = _mm_sub_pd(iy1,jy0);
 451             dz10             = _mm_sub_pd(iz1,jz0);
 452             dx20             = _mm_sub_pd(ix2,jx0);
 453             dy20             = _mm_sub_pd(iy2,jy0);
 454             dz20             = _mm_sub_pd(iz2,jz0);
 455
 456             /* Calculate squared distance and things based on it */
 457             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
 458             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
 459             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
 460
 461             rinv00           = gmx_mm_invsqrt_pd(rsq00);
 462             rinv10           = gmx_mm_invsqrt_pd(rsq10);
 463             rinv20           = gmx_mm_invsqrt_pd(rsq20);
 464
 465             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
 466             rinvsq10         = _mm_mul_pd(rinv10,rinv10);
 467             rinvsq20         = _mm_mul_pd(rinv20,rinv20);
 468
 469             /* Load parameters for j particles */
 470             jq0              = _mm_load_sd(charge+jnrA+0);
 471             vdwjidx0A        = 2*vdwtype[jnrA+0];
 472
 473             fjx0             = _mm_setzero_pd();
 474             fjy0             = _mm_setzero_pd();
 475             fjz0             = _mm_setzero_pd();
 476
 477             /**************************
 478              * CALCULATE INTERACTIONS *
 479              **************************/
 480
 481             if (gmx_mm_any_lt(rsq00,rcutoff2))
 482             {
 483
 484             r00              = _mm_mul_pd(rsq00,rinv00);
 485
 486             /* Compute parameters for interactions between i and j atoms */
 487             qq00             = _mm_mul_pd(iq0,jq0);
 488             gmx_mm_load_1pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,&c6_00,&c12_00);
 489             c6grid_00       = gmx_mm_load_1real_pd(vdwgridparam+vdwioffset0+vdwjidx0A);
 490
 491             /* EWALD ELECTROSTATICS */
 492
 493             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 494             ewrt             = _mm_mul_pd(r00,ewtabscale);
 495             ewitab           = _mm_cvttpd_epi32(ewrt);
 496 #ifdef __XOP__
 497             eweps            = _mm_frcz_pd(ewrt);
 498 #else
 499             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 500 #endif
 501             twoeweps         = _mm_add_pd(eweps,eweps);
 502             ewitab           = _mm_slli_epi32(ewitab,2);
 503             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 504             ewtabD           = _mm_setzero_pd();
 505             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 506             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 507             ewtabFn          = _mm_setzero_pd();
 508             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 509             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 510             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 511             velec            = _mm_mul_pd(qq00,_mm_sub_pd(_mm_sub_pd(rinv00,sh_ewald),velec));
 512             felec            = _mm_mul_pd(_mm_mul_pd(qq00,rinv00),_mm_sub_pd(rinvsq00,felec));
 513
 514             /* Analytical LJ-PME */
 515             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
 516             ewcljrsq         = _mm_mul_pd(ewclj2,rsq00);
 517             ewclj6           = _mm_mul_pd(ewclj2,_mm_mul_pd(ewclj2,ewclj2));
 518             exponent         = gmx_simd_exp_d(ewcljrsq);
 519             /* poly = exp(-(beta*r)^2) * (1 + (beta*r)^2 + (beta*r)^4 /2) */
 520             poly             = _mm_mul_pd(exponent,_mm_macc_pd(_mm_mul_pd(ewcljrsq,ewcljrsq),one_half,_mm_sub_pd(one,ewcljrsq)));
 521             /* vvdw6 = [C6 - C6grid * (1-poly)]/r6 */
 522             vvdw6            = _mm_mul_pd(_mm_macc_pd(-c6grid_00,_mm_sub_pd(one,poly),c6_00),rinvsix);
 523             vvdw12           = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
 524             vvdw             = _mm_msub_pd(_mm_nmacc_pd(c12_00,_mm_mul_pd(sh_vdw_invrcut6,sh_vdw_invrcut6),vvdw12),one_twelfth,
 525                                _mm_mul_pd(_mm_sub_pd(vvdw6,_mm_macc_pd(c6grid_00,sh_lj_ewald,_mm_mul_pd(c6_00,sh_vdw_invrcut6))),one_sixth));
 526             /* fvdw = vvdw12/r - (vvdw6/r + (C6grid * exponent * beta^6)/r) */
 527             fvdw             = _mm_mul_pd(_mm_add_pd(vvdw12,_mm_msub_pd(_mm_mul_pd(c6grid_00,one_sixth),_mm_mul_pd(exponent,ewclj6),vvdw6)),rinvsq00);
 528
 529             cutoff_mask      = _mm_cmplt_pd(rsq00,rcutoff2);
 530
 531             /* Update potential sum for this i atom from the interaction with this j atom. */
 532             velec            = _mm_and_pd(velec,cutoff_mask);
 533             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
 534             velecsum         = _mm_add_pd(velecsum,velec);
 535             vvdw             = _mm_and_pd(vvdw,cutoff_mask);
 536             vvdw             = _mm_unpacklo_pd(vvdw,_mm_setzero_pd());
 537             vvdwsum          = _mm_add_pd(vvdwsum,vvdw);
 538
 539             fscal            = _mm_add_pd(felec,fvdw);
 540
 541             fscal            = _mm_and_pd(fscal,cutoff_mask);
 542
 543             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
 544
 545             /* Update vectorial force */
 546             fix0             = _mm_macc_pd(dx00,fscal,fix0);
 547             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
 548             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
 549
 550             fjx0             = _mm_macc_pd(dx00,fscal,fjx0);
 551             fjy0             = _mm_macc_pd(dy00,fscal,fjy0);
 552             fjz0             = _mm_macc_pd(dz00,fscal,fjz0);
 553
 554             }
 555
 556             /**************************
 557              * CALCULATE INTERACTIONS *
 558              **************************/
 559
 560             if (gmx_mm_any_lt(rsq10,rcutoff2))
 561             {
 562
 563             r10              = _mm_mul_pd(rsq10,rinv10);
 564
 565             /* Compute parameters for interactions between i and j atoms */
 566             qq10             = _mm_mul_pd(iq1,jq0);
 567
 568             /* EWALD ELECTROSTATICS */
 569
 570             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 571             ewrt             = _mm_mul_pd(r10,ewtabscale);
 572             ewitab           = _mm_cvttpd_epi32(ewrt);
 573 #ifdef __XOP__
 574             eweps            = _mm_frcz_pd(ewrt);
 575 #else
 576             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 577 #endif
 578             twoeweps         = _mm_add_pd(eweps,eweps);
 579             ewitab           = _mm_slli_epi32(ewitab,2);
 580             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 581             ewtabD           = _mm_setzero_pd();
 582             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 583             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 584             ewtabFn          = _mm_setzero_pd();
 585             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 586             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 587             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 588             velec            = _mm_mul_pd(qq10,_mm_sub_pd(_mm_sub_pd(rinv10,sh_ewald),velec));
 589             felec            = _mm_mul_pd(_mm_mul_pd(qq10,rinv10),_mm_sub_pd(rinvsq10,felec));
 590
 591             cutoff_mask      = _mm_cmplt_pd(rsq10,rcutoff2);
 592
 593             /* Update potential sum for this i atom from the interaction with this j atom. */
 594             velec            = _mm_and_pd(velec,cutoff_mask);
 595             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
 596             velecsum         = _mm_add_pd(velecsum,velec);
 597
 598             fscal            = felec;
 599
 600             fscal            = _mm_and_pd(fscal,cutoff_mask);
 601
 602             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
 603
 604             /* Update vectorial force */
 605             fix1             = _mm_macc_pd(dx10,fscal,fix1);
 606             fiy1             = _mm_macc_pd(dy10,fscal,fiy1);
 607             fiz1             = _mm_macc_pd(dz10,fscal,fiz1);
 608
 609             fjx0             = _mm_macc_pd(dx10,fscal,fjx0);
 610             fjy0             = _mm_macc_pd(dy10,fscal,fjy0);
 611             fjz0             = _mm_macc_pd(dz10,fscal,fjz0);
 612
 613             }
 614
 615             /**************************
 616              * CALCULATE INTERACTIONS *
 617              **************************/
 618
 619             if (gmx_mm_any_lt(rsq20,rcutoff2))
 620             {
 621
 622             r20              = _mm_mul_pd(rsq20,rinv20);
 623
 624             /* Compute parameters for interactions between i and j atoms */
 625             qq20             = _mm_mul_pd(iq2,jq0);
 626
 627             /* EWALD ELECTROSTATICS */
 628
 629             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 630             ewrt             = _mm_mul_pd(r20,ewtabscale);
 631             ewitab           = _mm_cvttpd_epi32(ewrt);
 632 #ifdef __XOP__
 633             eweps            = _mm_frcz_pd(ewrt);
 634 #else
 635             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 636 #endif
 637             twoeweps         = _mm_add_pd(eweps,eweps);
 638             ewitab           = _mm_slli_epi32(ewitab,2);
 639             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 640             ewtabD           = _mm_setzero_pd();
 641             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 642             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 643             ewtabFn          = _mm_setzero_pd();
 644             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 645             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 646             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 647             velec            = _mm_mul_pd(qq20,_mm_sub_pd(_mm_sub_pd(rinv20,sh_ewald),velec));
 648             felec            = _mm_mul_pd(_mm_mul_pd(qq20,rinv20),_mm_sub_pd(rinvsq20,felec));
 649
 650             cutoff_mask      = _mm_cmplt_pd(rsq20,rcutoff2);
 651
 652             /* Update potential sum for this i atom from the interaction with this j atom. */
 653             velec            = _mm_and_pd(velec,cutoff_mask);
 654             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
 655             velecsum         = _mm_add_pd(velecsum,velec);
 656
 657             fscal            = felec;
 658
 659             fscal            = _mm_and_pd(fscal,cutoff_mask);
 660
 661             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
 662
 663             /* Update vectorial force */
 664             fix2             = _mm_macc_pd(dx20,fscal,fix2);
 665             fiy2             = _mm_macc_pd(dy20,fscal,fiy2);
 666             fiz2             = _mm_macc_pd(dz20,fscal,fiz2);
 667
 668             fjx0             = _mm_macc_pd(dx20,fscal,fjx0);
 669             fjy0             = _mm_macc_pd(dy20,fscal,fjy0);
 670             fjz0             = _mm_macc_pd(dz20,fscal,fjz0);
 671
 672             }
 673
 674             gmx_mm_decrement_1rvec_1ptr_swizzle_pd(f+j_coord_offsetA,fjx0,fjy0,fjz0);
 675
 676             /* Inner loop uses 179 flops */
 677         }
 678
 679         /* End of innermost loop */
 680
 681         gmx_mm_update_iforce_3atom_swizzle_pd(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,
 682                                               f+i_coord_offset,fshift+i_shift_offset);
 683
 684         ggid                        = gid[iidx];
 685         /* Update potential energies */
 686         gmx_mm_update_1pot_pd(velecsum,kernel_data->energygrp_elec+ggid);
 687         gmx_mm_update_1pot_pd(vvdwsum,kernel_data->energygrp_vdw+ggid);
 688
 689         /* Increment number of inner iterations */
 690         inneriter                  += j_index_end - j_index_start;
 691
 692         /* Outer loop uses 20 flops */
 693     }
 694
 695     /* Increment number of outer iterations */
 696     outeriter        += nri;
 697
 698     /* Update outer/inner flops */
 699
 700     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W3_VF,outeriter*20 + inneriter*179);
 701 }
 702 /*
 703  * Gromacs nonbonded kernel:   nb_kernel_ElecEwSh_VdwLJEwSh_GeomW3P1_F_avx_128_fma_double
 704  * Electrostatics interaction: Ewald
 705  * VdW interaction:            LJEwald
 706  * Geometry:                   Water3-Particle
 707  * Calculate force/pot:        Force
 708  */
 709 void
 710 nb_kernel_ElecEwSh_VdwLJEwSh_GeomW3P1_F_avx_128_fma_double
 711                     (t_nblist                    * gmx_restrict       nlist,
 712                      rvec                        * gmx_restrict          xx,
 713                      rvec                        * gmx_restrict          ff,
 714                      t_forcerec                  * gmx_restrict          fr,
 715                      t_mdatoms                   * gmx_restrict     mdatoms,
 716                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
 717                      t_nrnb                      * gmx_restrict        nrnb)
 718 {
 719     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
 720      * just 0 for non-waters.
 721      * Suffixes A,B refer to j loop unrolling done with SSE double precision, e.g. for the two different
 722      * jnr indices corresponding to data put in the four positions in the SIMD register.
 723      */
 724     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
 725     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
 726     int              jnrA,jnrB;
 727     int              j_coord_offsetA,j_coord_offsetB;
 728     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
 729     real             rcutoff_scalar;
 730     real             *shiftvec,*fshift,*x,*f;
 731     __m128d          tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
 732     int              vdwioffset0;
 733     __m128d          ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
 734     int              vdwioffset1;
 735     __m128d          ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
 736     int              vdwioffset2;
 737     __m128d          ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
 738     int              vdwjidx0A,vdwjidx0B;
 739     __m128d          jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
 740     __m128d          dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
 741     __m128d          dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
 742     __m128d          dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
 743     __m128d          velec,felec,velecsum,facel,crf,krf,krf2;
 744     real             *charge;
 745     int              nvdwtype;
 746     __m128d          rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
 747     int              *vdwtype;
 748     real             *vdwparam;
 749     __m128d          one_sixth   = _mm_set1_pd(1.0/6.0);
 750     __m128d          one_twelfth = _mm_set1_pd(1.0/12.0);
 751     __m128d           c6grid_00;
 752     __m128d           c6grid_10;
 753     __m128d           c6grid_20;
 754     real             *vdwgridparam;
 755     __m128d           ewclj,ewclj2,ewclj6,ewcljrsq,poly,exponent,f6A,f6B,sh_lj_ewald;
 756     __m128d           one_half  = _mm_set1_pd(0.5);
 757     __m128d           minus_one = _mm_set1_pd(-1.0);
 758     __m128i          ewitab;
 759     __m128d          ewtabscale,eweps,twoeweps,sh_ewald,ewrt,ewtabhalfspace,ewtabF,ewtabFn,ewtabD,ewtabV;
 760     real             *ewtab;
 761     __m128d          dummy_mask,cutoff_mask;
 762     __m128d          signbit   = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
 763     __m128d          one     = _mm_set1_pd(1.0);
 764     __m128d          two     = _mm_set1_pd(2.0);
 765     x                = xx[0];
 766     f                = ff[0];
 767
 768     nri              = nlist->nri;
 769     iinr             = nlist->iinr;
 770     jindex           = nlist->jindex;
 771     jjnr             = nlist->jjnr;
 772     shiftidx         = nlist->shift;
 773     gid              = nlist->gid;
 774     shiftvec         = fr->shift_vec[0];
 775     fshift           = fr->fshift[0];
 776     facel            = _mm_set1_pd(fr->epsfac);
 777     charge           = mdatoms->chargeA;
 778     nvdwtype         = fr->ntype;
 779     vdwparam         = fr->nbfp;
 780     vdwtype          = mdatoms->typeA;
 781     vdwgridparam     = fr->ljpme_c6grid;
 782     sh_lj_ewald      = _mm_set1_pd(fr->ic->sh_lj_ewald);
 783     ewclj            = _mm_set1_pd(fr->ewaldcoeff_lj);
 784     ewclj2           = _mm_mul_pd(minus_one,_mm_mul_pd(ewclj,ewclj));
 785
 786     sh_ewald         = _mm_set1_pd(fr->ic->sh_ewald);
 787     ewtab            = fr->ic->tabq_coul_F;
 788     ewtabscale       = _mm_set1_pd(fr->ic->tabq_scale);
 789     ewtabhalfspace   = _mm_set1_pd(0.5/fr->ic->tabq_scale);
 790
 791     /* Setup water-specific parameters */
 792     inr              = nlist->iinr[0];
 793     iq0              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+0]));
 794     iq1              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+1]));
 795     iq2              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+2]));
 796     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
 797
 798     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
 799     rcutoff_scalar   = fr->rcoulomb;
 800     rcutoff          = _mm_set1_pd(rcutoff_scalar);
 801     rcutoff2         = _mm_mul_pd(rcutoff,rcutoff);
 802
 803     sh_vdw_invrcut6  = _mm_set1_pd(fr->ic->sh_invrc6);
 804     rvdw             = _mm_set1_pd(fr->rvdw);
 805
 806     /* Avoid stupid compiler warnings */
 807     jnrA = jnrB = 0;
 808     j_coord_offsetA = 0;
 809     j_coord_offsetB = 0;
 810
 811     outeriter        = 0;
 812     inneriter        = 0;
 813
 814     /* Start outer loop over neighborlists */
 815     for(iidx=0; iidx<nri; iidx++)
 816     {
 817         /* Load shift vector for this list */
 818         i_shift_offset   = DIM*shiftidx[iidx];
 819
 820         /* Load limits for loop over neighbors */
 821         j_index_start    = jindex[iidx];
 822         j_index_end      = jindex[iidx+1];
 823
 824         /* Get outer coordinate index */
 825         inr              = iinr[iidx];
 826         i_coord_offset   = DIM*inr;
 827
 828         /* Load i particle coords and add shift vector */
 829         gmx_mm_load_shift_and_3rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,
 830                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2);
 831
 832         fix0             = _mm_setzero_pd();
 833         fiy0             = _mm_setzero_pd();
 834         fiz0             = _mm_setzero_pd();
 835         fix1             = _mm_setzero_pd();
 836         fiy1             = _mm_setzero_pd();
 837         fiz1             = _mm_setzero_pd();
 838         fix2             = _mm_setzero_pd();
 839         fiy2             = _mm_setzero_pd();
 840         fiz2             = _mm_setzero_pd();
 841
 842         /* Start inner kernel loop */
 843         for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
 844         {
 845
 846             /* Get j neighbor index, and coordinate index */
 847             jnrA             = jjnr[jidx];
 848             jnrB             = jjnr[jidx+1];
 849             j_coord_offsetA  = DIM*jnrA;
 850             j_coord_offsetB  = DIM*jnrB;
 851
 852             /* load j atom coordinates */
 853             gmx_mm_load_1rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
 854                                               &jx0,&jy0,&jz0);
 855
 856             /* Calculate displacement vector */
 857             dx00             = _mm_sub_pd(ix0,jx0);
 858             dy00             = _mm_sub_pd(iy0,jy0);
 859             dz00             = _mm_sub_pd(iz0,jz0);
 860             dx10             = _mm_sub_pd(ix1,jx0);
 861             dy10             = _mm_sub_pd(iy1,jy0);
 862             dz10             = _mm_sub_pd(iz1,jz0);
 863             dx20             = _mm_sub_pd(ix2,jx0);
 864             dy20             = _mm_sub_pd(iy2,jy0);
 865             dz20             = _mm_sub_pd(iz2,jz0);
 866
 867             /* Calculate squared distance and things based on it */
 868             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
 869             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
 870             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
 871
 872             rinv00           = gmx_mm_invsqrt_pd(rsq00);
 873             rinv10           = gmx_mm_invsqrt_pd(rsq10);
 874             rinv20           = gmx_mm_invsqrt_pd(rsq20);
 875
 876             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
 877             rinvsq10         = _mm_mul_pd(rinv10,rinv10);
 878             rinvsq20         = _mm_mul_pd(rinv20,rinv20);
 879
 880             /* Load parameters for j particles */
 881             jq0              = gmx_mm_load_2real_swizzle_pd(charge+jnrA+0,charge+jnrB+0);
 882             vdwjidx0A        = 2*vdwtype[jnrA+0];
 883             vdwjidx0B        = 2*vdwtype[jnrB+0];
 884
 885             fjx0             = _mm_setzero_pd();
 886             fjy0             = _mm_setzero_pd();
 887             fjz0             = _mm_setzero_pd();
 888
 889             /**************************
 890              * CALCULATE INTERACTIONS *
 891              **************************/
 892
 893             if (gmx_mm_any_lt(rsq00,rcutoff2))
 894             {
 895
 896             r00              = _mm_mul_pd(rsq00,rinv00);
 897
 898             /* Compute parameters for interactions between i and j atoms */
 899             qq00             = _mm_mul_pd(iq0,jq0);
 900             gmx_mm_load_2pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,
 901                                          vdwparam+vdwioffset0+vdwjidx0B,&c6_00,&c12_00);
 902             c6grid_00       = gmx_mm_load_2real_swizzle_pd(vdwgridparam+vdwioffset0+vdwjidx0A,
 903                                                                vdwgridparam+vdwioffset0+vdwjidx0B);
 904
 905             /* EWALD ELECTROSTATICS */
 906
 907             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 908             ewrt             = _mm_mul_pd(r00,ewtabscale);
 909             ewitab           = _mm_cvttpd_epi32(ewrt);
 910 #ifdef __XOP__
 911             eweps            = _mm_frcz_pd(ewrt);
 912 #else
 913             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 914 #endif
 915             twoeweps         = _mm_add_pd(eweps,eweps);
 916             gmx_mm_load_2pair_swizzle_pd(ewtab+_mm_extract_epi32(ewitab,0),ewtab+_mm_extract_epi32(ewitab,1),
 917                                          &ewtabF,&ewtabFn);
 918             felec            = _mm_macc_pd(eweps,ewtabFn,_mm_mul_pd( _mm_sub_pd(one,eweps),ewtabF));
 919             felec            = _mm_mul_pd(_mm_mul_pd(qq00,rinv00),_mm_sub_pd(rinvsq00,felec));
 920
 921             /* Analytical LJ-PME */
 922             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
 923             ewcljrsq         = _mm_mul_pd(ewclj2,rsq00);
 924             ewclj6           = _mm_mul_pd(ewclj2,_mm_mul_pd(ewclj2,ewclj2));
 925             exponent         = gmx_simd_exp_d(ewcljrsq);
 926             /* poly = exp(-(beta*r)^2) * (1 + (beta*r)^2 + (beta*r)^4 /2) */
 927             poly             = _mm_mul_pd(exponent,_mm_macc_pd(_mm_mul_pd(ewcljrsq,ewcljrsq),one_half,_mm_sub_pd(one,ewcljrsq)));
 928             /* f6A = 6 * C6grid * (1 - poly) */
 929             f6A              = _mm_mul_pd(c6grid_00,_mm_sub_pd(one,poly));
 930             /* f6B = C6grid * exponent * beta^6 */
 931             f6B              = _mm_mul_pd(_mm_mul_pd(c6grid_00,one_sixth),_mm_mul_pd(exponent,ewclj6));
 932             /* fvdw = 12*C12/r13 - ((6*C6 - f6A)/r6 + f6B)/r */
 933             fvdw              = _mm_mul_pd(_mm_macc_pd(_mm_msub_pd(c12_00,rinvsix,_mm_sub_pd(c6_00,f6A)),rinvsix,f6B),rinvsq00);
 934
 935             cutoff_mask      = _mm_cmplt_pd(rsq00,rcutoff2);
 936
 937             fscal            = _mm_add_pd(felec,fvdw);
 938
 939             fscal            = _mm_and_pd(fscal,cutoff_mask);
 940
 941             /* Update vectorial force */
 942             fix0             = _mm_macc_pd(dx00,fscal,fix0);
 943             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
 944             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
 945
 946             fjx0             = _mm_macc_pd(dx00,fscal,fjx0);
 947             fjy0             = _mm_macc_pd(dy00,fscal,fjy0);
 948             fjz0             = _mm_macc_pd(dz00,fscal,fjz0);
 949
 950             }
 951
 952             /**************************
 953              * CALCULATE INTERACTIONS *
 954              **************************/
 955
 956             if (gmx_mm_any_lt(rsq10,rcutoff2))
 957             {
 958
 959             r10              = _mm_mul_pd(rsq10,rinv10);
 960
 961             /* Compute parameters for interactions between i and j atoms */
 962             qq10             = _mm_mul_pd(iq1,jq0);
 963
 964             /* EWALD ELECTROSTATICS */
 965
 966             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 967             ewrt             = _mm_mul_pd(r10,ewtabscale);
 968             ewitab           = _mm_cvttpd_epi32(ewrt);
 969 #ifdef __XOP__
 970             eweps            = _mm_frcz_pd(ewrt);
 971 #else
 972             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 973 #endif
 974             twoeweps         = _mm_add_pd(eweps,eweps);
 975             gmx_mm_load_2pair_swizzle_pd(ewtab+_mm_extract_epi32(ewitab,0),ewtab+_mm_extract_epi32(ewitab,1),
 976                                          &ewtabF,&ewtabFn);
 977             felec            = _mm_macc_pd(eweps,ewtabFn,_mm_mul_pd( _mm_sub_pd(one,eweps),ewtabF));
 978             felec            = _mm_mul_pd(_mm_mul_pd(qq10,rinv10),_mm_sub_pd(rinvsq10,felec));
 979
 980             cutoff_mask      = _mm_cmplt_pd(rsq10,rcutoff2);
 981
 982             fscal            = felec;
 983
 984             fscal            = _mm_and_pd(fscal,cutoff_mask);
 985
 986             /* Update vectorial force */
 987             fix1             = _mm_macc_pd(dx10,fscal,fix1);
 988             fiy1             = _mm_macc_pd(dy10,fscal,fiy1);
 989             fiz1             = _mm_macc_pd(dz10,fscal,fiz1);
 990
 991             fjx0             = _mm_macc_pd(dx10,fscal,fjx0);
 992             fjy0             = _mm_macc_pd(dy10,fscal,fjy0);
 993             fjz0             = _mm_macc_pd(dz10,fscal,fjz0);
 994
 995             }
 996
 997             /**************************
 998              * CALCULATE INTERACTIONS *
 999              **************************/
1000
1001             if (gmx_mm_any_lt(rsq20,rcutoff2))
1002             {
1003
1004             r20              = _mm_mul_pd(rsq20,rinv20);
1005
1006             /* Compute parameters for interactions between i and j atoms */
1007             qq20             = _mm_mul_pd(iq2,jq0);
1008
1009             /* EWALD ELECTROSTATICS */
1010
1011             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1012             ewrt             = _mm_mul_pd(r20,ewtabscale);
1013             ewitab           = _mm_cvttpd_epi32(ewrt);
1014 #ifdef __XOP__
1015             eweps            = _mm_frcz_pd(ewrt);
1016 #else
1017             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1018 #endif
1019             twoeweps         = _mm_add_pd(eweps,eweps);
1020             gmx_mm_load_2pair_swizzle_pd(ewtab+_mm_extract_epi32(ewitab,0),ewtab+_mm_extract_epi32(ewitab,1),
1021                                          &ewtabF,&ewtabFn);
1022             felec            = _mm_macc_pd(eweps,ewtabFn,_mm_mul_pd( _mm_sub_pd(one,eweps),ewtabF));
1023             felec            = _mm_mul_pd(_mm_mul_pd(qq20,rinv20),_mm_sub_pd(rinvsq20,felec));
1024
1025             cutoff_mask      = _mm_cmplt_pd(rsq20,rcutoff2);
1026
1027             fscal            = felec;
1028
1029             fscal            = _mm_and_pd(fscal,cutoff_mask);
1030
1031             /* Update vectorial force */
1032             fix2             = _mm_macc_pd(dx20,fscal,fix2);
1033             fiy2             = _mm_macc_pd(dy20,fscal,fiy2);
1034             fiz2             = _mm_macc_pd(dz20,fscal,fiz2);
1035
1036             fjx0             = _mm_macc_pd(dx20,fscal,fjx0);
1037             fjy0             = _mm_macc_pd(dy20,fscal,fjy0);
1038             fjz0             = _mm_macc_pd(dz20,fscal,fjz0);
1039
1040             }
1041
1042             gmx_mm_decrement_1rvec_2ptr_swizzle_pd(f+j_coord_offsetA,f+j_coord_offsetB,fjx0,fjy0,fjz0);
1043
1044             /* Inner loop uses 150 flops */
1045         }
1046
1047         if(jidx<j_index_end)
1048         {
1049
1050             jnrA             = jjnr[jidx];
1051             j_coord_offsetA  = DIM*jnrA;
1052
1053             /* load j atom coordinates */
1054             gmx_mm_load_1rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
1055                                               &jx0,&jy0,&jz0);
1056
1057             /* Calculate displacement vector */
1058             dx00             = _mm_sub_pd(ix0,jx0);
1059             dy00             = _mm_sub_pd(iy0,jy0);
1060             dz00             = _mm_sub_pd(iz0,jz0);
1061             dx10             = _mm_sub_pd(ix1,jx0);
1062             dy10             = _mm_sub_pd(iy1,jy0);
1063             dz10             = _mm_sub_pd(iz1,jz0);
1064             dx20             = _mm_sub_pd(ix2,jx0);
1065             dy20             = _mm_sub_pd(iy2,jy0);
1066             dz20             = _mm_sub_pd(iz2,jz0);
1067
1068             /* Calculate squared distance and things based on it */
1069             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
1070             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
1071             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
1072
1073             rinv00           = gmx_mm_invsqrt_pd(rsq00);
1074             rinv10           = gmx_mm_invsqrt_pd(rsq10);
1075             rinv20           = gmx_mm_invsqrt_pd(rsq20);
1076
1077             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
1078             rinvsq10         = _mm_mul_pd(rinv10,rinv10);
1079             rinvsq20         = _mm_mul_pd(rinv20,rinv20);
1080
1081             /* Load parameters for j particles */
1082             jq0              = _mm_load_sd(charge+jnrA+0);
1083             vdwjidx0A        = 2*vdwtype[jnrA+0];
1084
1085             fjx0             = _mm_setzero_pd();
1086             fjy0             = _mm_setzero_pd();
1087             fjz0             = _mm_setzero_pd();
1088
1089             /**************************
1090              * CALCULATE INTERACTIONS *
1091              **************************/
1092
1093             if (gmx_mm_any_lt(rsq00,rcutoff2))
1094             {
1095
1096             r00              = _mm_mul_pd(rsq00,rinv00);
1097
1098             /* Compute parameters for interactions between i and j atoms */
1099             qq00             = _mm_mul_pd(iq0,jq0);
1100             gmx_mm_load_1pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,&c6_00,&c12_00);
1101             c6grid_00       = gmx_mm_load_1real_pd(vdwgridparam+vdwioffset0+vdwjidx0A);
1102
1103             /* EWALD ELECTROSTATICS */
1104
1105             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1106             ewrt             = _mm_mul_pd(r00,ewtabscale);
1107             ewitab           = _mm_cvttpd_epi32(ewrt);
1108 #ifdef __XOP__
1109             eweps            = _mm_frcz_pd(ewrt);
1110 #else
1111             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1112 #endif
1113             twoeweps         = _mm_add_pd(eweps,eweps);
1114             gmx_mm_load_1pair_swizzle_pd(ewtab+_mm_extract_epi32(ewitab,0),&ewtabF,&ewtabFn);
1115             felec            = _mm_macc_pd(eweps,ewtabFn,_mm_mul_pd( _mm_sub_pd(one,eweps),ewtabF));
1116             felec            = _mm_mul_pd(_mm_mul_pd(qq00,rinv00),_mm_sub_pd(rinvsq00,felec));
1117
1118             /* Analytical LJ-PME */
1119             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
1120             ewcljrsq         = _mm_mul_pd(ewclj2,rsq00);
1121             ewclj6           = _mm_mul_pd(ewclj2,_mm_mul_pd(ewclj2,ewclj2));
1122             exponent         = gmx_simd_exp_d(ewcljrsq);
1123             /* poly = exp(-(beta*r)^2) * (1 + (beta*r)^2 + (beta*r)^4 /2) */
1124             poly             = _mm_mul_pd(exponent,_mm_macc_pd(_mm_mul_pd(ewcljrsq,ewcljrsq),one_half,_mm_sub_pd(one,ewcljrsq)));
1125             /* f6A = 6 * C6grid * (1 - poly) */
1126             f6A              = _mm_mul_pd(c6grid_00,_mm_sub_pd(one,poly));
1127             /* f6B = C6grid * exponent * beta^6 */
1128             f6B              = _mm_mul_pd(_mm_mul_pd(c6grid_00,one_sixth),_mm_mul_pd(exponent,ewclj6));
1129             /* fvdw = 12*C12/r13 - ((6*C6 - f6A)/r6 + f6B)/r */
1130             fvdw              = _mm_mul_pd(_mm_macc_pd(_mm_msub_pd(c12_00,rinvsix,_mm_sub_pd(c6_00,f6A)),rinvsix,f6B),rinvsq00);
1131
1132             cutoff_mask      = _mm_cmplt_pd(rsq00,rcutoff2);
1133
1134             fscal            = _mm_add_pd(felec,fvdw);
1135
1136             fscal            = _mm_and_pd(fscal,cutoff_mask);
1137
1138             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1139
1140             /* Update vectorial force */
1141             fix0             = _mm_macc_pd(dx00,fscal,fix0);
1142             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
1143             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
1144
1145             fjx0             = _mm_macc_pd(dx00,fscal,fjx0);
1146             fjy0             = _mm_macc_pd(dy00,fscal,fjy0);
1147             fjz0             = _mm_macc_pd(dz00,fscal,fjz0);
1148
1149             }
1150
1151             /**************************
1152              * CALCULATE INTERACTIONS *
1153              **************************/
1154
1155             if (gmx_mm_any_lt(rsq10,rcutoff2))
1156             {
1157
1158             r10              = _mm_mul_pd(rsq10,rinv10);
1159
1160             /* Compute parameters for interactions between i and j atoms */
1161             qq10             = _mm_mul_pd(iq1,jq0);
1162
1163             /* EWALD ELECTROSTATICS */
1164
1165             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1166             ewrt             = _mm_mul_pd(r10,ewtabscale);
1167             ewitab           = _mm_cvttpd_epi32(ewrt);
1168 #ifdef __XOP__
1169             eweps            = _mm_frcz_pd(ewrt);
1170 #else
1171             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1172 #endif
1173             twoeweps         = _mm_add_pd(eweps,eweps);
1174             gmx_mm_load_1pair_swizzle_pd(ewtab+_mm_extract_epi32(ewitab,0),&ewtabF,&ewtabFn);
1175             felec            = _mm_macc_pd(eweps,ewtabFn,_mm_mul_pd( _mm_sub_pd(one,eweps),ewtabF));
1176             felec            = _mm_mul_pd(_mm_mul_pd(qq10,rinv10),_mm_sub_pd(rinvsq10,felec));
1177
1178             cutoff_mask      = _mm_cmplt_pd(rsq10,rcutoff2);
1179
1180             fscal            = felec;
1181
1182             fscal            = _mm_and_pd(fscal,cutoff_mask);
1183
1184             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1185
1186             /* Update vectorial force */
1187             fix1             = _mm_macc_pd(dx10,fscal,fix1);
1188             fiy1             = _mm_macc_pd(dy10,fscal,fiy1);
1189             fiz1             = _mm_macc_pd(dz10,fscal,fiz1);
1190
1191             fjx0             = _mm_macc_pd(dx10,fscal,fjx0);
1192             fjy0             = _mm_macc_pd(dy10,fscal,fjy0);
1193             fjz0             = _mm_macc_pd(dz10,fscal,fjz0);
1194
1195             }
1196
1197             /**************************
1198              * CALCULATE INTERACTIONS *
1199              **************************/
1200
1201             if (gmx_mm_any_lt(rsq20,rcutoff2))
1202             {
1203
1204             r20              = _mm_mul_pd(rsq20,rinv20);
1205
1206             /* Compute parameters for interactions between i and j atoms */
1207             qq20             = _mm_mul_pd(iq2,jq0);
1208
1209             /* EWALD ELECTROSTATICS */
1210
1211             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1212             ewrt             = _mm_mul_pd(r20,ewtabscale);
1213             ewitab           = _mm_cvttpd_epi32(ewrt);
1214 #ifdef __XOP__
1215             eweps            = _mm_frcz_pd(ewrt);
1216 #else
1217             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1218 #endif
1219             twoeweps         = _mm_add_pd(eweps,eweps);
1220             gmx_mm_load_1pair_swizzle_pd(ewtab+_mm_extract_epi32(ewitab,0),&ewtabF,&ewtabFn);
1221             felec            = _mm_macc_pd(eweps,ewtabFn,_mm_mul_pd( _mm_sub_pd(one,eweps),ewtabF));
1222             felec            = _mm_mul_pd(_mm_mul_pd(qq20,rinv20),_mm_sub_pd(rinvsq20,felec));
1223
1224             cutoff_mask      = _mm_cmplt_pd(rsq20,rcutoff2);
1225
1226             fscal            = felec;
1227
1228             fscal            = _mm_and_pd(fscal,cutoff_mask);
1229
1230             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1231
1232             /* Update vectorial force */
1233             fix2             = _mm_macc_pd(dx20,fscal,fix2);
1234             fiy2             = _mm_macc_pd(dy20,fscal,fiy2);
1235             fiz2             = _mm_macc_pd(dz20,fscal,fiz2);
1236
1237             fjx0             = _mm_macc_pd(dx20,fscal,fjx0);
1238             fjy0             = _mm_macc_pd(dy20,fscal,fjy0);
1239             fjz0             = _mm_macc_pd(dz20,fscal,fjz0);
1240
1241             }
1242
1243             gmx_mm_decrement_1rvec_1ptr_swizzle_pd(f+j_coord_offsetA,fjx0,fjy0,fjz0);
1244
1245             /* Inner loop uses 150 flops */
1246         }
1247
1248         /* End of innermost loop */
1249
1250         gmx_mm_update_iforce_3atom_swizzle_pd(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,
1251                                               f+i_coord_offset,fshift+i_shift_offset);
1252
1253         /* Increment number of inner iterations */
1254         inneriter                  += j_index_end - j_index_start;
1255
1256         /* Outer loop uses 18 flops */
1257     }
1258
1259     /* Increment number of outer iterations */
1260     outeriter        += nri;
1261
1262     /* Update outer/inner flops */
1263
1264     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W3_F,outeriter*18 + inneriter*150);
1265 }