src/gromacs/gmxlib/nonbonded/nb_kernel_avx_128_fma_single/nb_kernel_ElecEwSh_VdwLJEwSh_GeomW3W3_avx_128_fma_single.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014,2015,2017, 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_single kernel generator.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "config.h"
  41
  42 #include <math.h>
  43
  44 #include "../nb_kernel.h"
  45 #include "gromacs/gmxlib/nrnb.h"
  46
  47 #include "kernelutil_x86_avx_128_fma_single.h"
  48
  49 /*
  50  * Gromacs nonbonded kernel:   nb_kernel_ElecEwSh_VdwLJEwSh_GeomW3W3_VF_avx_128_fma_single
  51  * Electrostatics interaction: Ewald
  52  * VdW interaction:            LJEwald
  53  * Geometry:                   Water3-Water3
  54  * Calculate force/pot:        PotentialAndForce
  55  */
  56 void
  57 nb_kernel_ElecEwSh_VdwLJEwSh_GeomW3W3_VF_avx_128_fma_single
  58                     (t_nblist                    * gmx_restrict       nlist,
  59                      rvec                        * gmx_restrict          xx,
  60                      rvec                        * gmx_restrict          ff,
  61                      struct t_forcerec           * gmx_restrict          fr,
  62                      t_mdatoms                   * gmx_restrict     mdatoms,
  63                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
  64                      t_nrnb                      * gmx_restrict        nrnb)
  65 {
  66     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
  67      * just 0 for non-waters.
  68      * Suffixes A,B,C,D refer to j loop unrolling done with AVX_128, e.g. for the four different
  69      * jnr indices corresponding to data put in the four positions in the SIMD register.
  70      */
  71     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
  72     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
  73     int              jnrA,jnrB,jnrC,jnrD;
  74     int              jnrlistA,jnrlistB,jnrlistC,jnrlistD;
  75     int              j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
  76     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
  77     real             rcutoff_scalar;
  78     real             *shiftvec,*fshift,*x,*f;
  79     real             *fjptrA,*fjptrB,*fjptrC,*fjptrD;
  80     real             scratch[4*DIM];
  81     __m128           fscal,rcutoff,rcutoff2,jidxall;
  82     int              vdwioffset0;
  83     __m128           ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
  84     int              vdwioffset1;
  85     __m128           ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
  86     int              vdwioffset2;
  87     __m128           ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
  88     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
  89     __m128           jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
  90     int              vdwjidx1A,vdwjidx1B,vdwjidx1C,vdwjidx1D;
  91     __m128           jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
  92     int              vdwjidx2A,vdwjidx2B,vdwjidx2C,vdwjidx2D;
  93     __m128           jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
  94     __m128           dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
  95     __m128           dx01,dy01,dz01,rsq01,rinv01,rinvsq01,r01,qq01,c6_01,c12_01;
  96     __m128           dx02,dy02,dz02,rsq02,rinv02,rinvsq02,r02,qq02,c6_02,c12_02;
  97     __m128           dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
  98     __m128           dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11;
  99     __m128           dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12;
 100     __m128           dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
 101     __m128           dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21;
 102     __m128           dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22;
 103     __m128           velec,felec,velecsum,facel,crf,krf,krf2;
 104     real             *charge;
 105     int              nvdwtype;
 106     __m128           rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
 107     int              *vdwtype;
 108     real             *vdwparam;
 109     __m128           one_sixth   = _mm_set1_ps(1.0/6.0);
 110     __m128           one_twelfth = _mm_set1_ps(1.0/12.0);
 111     __m128           c6grid_00;
 112     __m128           c6grid_01;
 113     __m128           c6grid_02;
 114     __m128           c6grid_10;
 115     __m128           c6grid_11;
 116     __m128           c6grid_12;
 117     __m128           c6grid_20;
 118     __m128           c6grid_21;
 119     __m128           c6grid_22;
 120     real             *vdwgridparam;
 121     __m128           ewclj,ewclj2,ewclj6,ewcljrsq,poly,exponent,f6A,f6B,sh_lj_ewald;
 122     __m128           one_half = _mm_set1_ps(0.5);
 123     __m128           minus_one = _mm_set1_ps(-1.0);
 124     __m128i          ewitab;
 125     __m128           ewtabscale,eweps,twoeweps,sh_ewald,ewrt,ewtabhalfspace,ewtabF,ewtabFn,ewtabD,ewtabV;
 126     __m128           beta,beta2,beta3,zeta2,pmecorrF,pmecorrV,rinv3;
 127     real             *ewtab;
 128     __m128           dummy_mask,cutoff_mask;
 129     __m128           signbit = _mm_castsi128_ps( _mm_set1_epi32(0x80000000) );
 130     __m128           one     = _mm_set1_ps(1.0);
 131     __m128           two     = _mm_set1_ps(2.0);
 132     x                = xx[0];
 133     f                = ff[0];
 134
 135     nri              = nlist->nri;
 136     iinr             = nlist->iinr;
 137     jindex           = nlist->jindex;
 138     jjnr             = nlist->jjnr;
 139     shiftidx         = nlist->shift;
 140     gid              = nlist->gid;
 141     shiftvec         = fr->shift_vec[0];
 142     fshift           = fr->fshift[0];
 143     facel            = _mm_set1_ps(fr->ic->epsfac);
 144     charge           = mdatoms->chargeA;
 145     nvdwtype         = fr->ntype;
 146     vdwparam         = fr->nbfp;
 147     vdwtype          = mdatoms->typeA;
 148     vdwgridparam     = fr->ljpme_c6grid;
 149     sh_lj_ewald      = _mm_set1_ps(fr->ic->sh_lj_ewald);
 150     ewclj            = _mm_set1_ps(fr->ic->ewaldcoeff_lj);
 151     ewclj2           = _mm_mul_ps(minus_one,_mm_mul_ps(ewclj,ewclj));
 152
 153     sh_ewald         = _mm_set1_ps(fr->ic->sh_ewald);
 154     beta             = _mm_set1_ps(fr->ic->ewaldcoeff_q);
 155     beta2            = _mm_mul_ps(beta,beta);
 156     beta3            = _mm_mul_ps(beta,beta2);
 157     ewtab            = fr->ic->tabq_coul_FDV0;
 158     ewtabscale       = _mm_set1_ps(fr->ic->tabq_scale);
 159     ewtabhalfspace   = _mm_set1_ps(0.5/fr->ic->tabq_scale);
 160
 161     /* Setup water-specific parameters */
 162     inr              = nlist->iinr[0];
 163     iq0              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+0]));
 164     iq1              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+1]));
 165     iq2              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+2]));
 166     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
 167
 168     jq0              = _mm_set1_ps(charge[inr+0]);
 169     jq1              = _mm_set1_ps(charge[inr+1]);
 170     jq2              = _mm_set1_ps(charge[inr+2]);
 171     vdwjidx0A        = 2*vdwtype[inr+0];
 172     qq00             = _mm_mul_ps(iq0,jq0);
 173     c6_00            = _mm_set1_ps(vdwparam[vdwioffset0+vdwjidx0A]);
 174     c12_00           = _mm_set1_ps(vdwparam[vdwioffset0+vdwjidx0A+1]);
 175     c6grid_00        = _mm_set1_ps(vdwgridparam[vdwioffset0+vdwjidx0A]);
 176     qq01             = _mm_mul_ps(iq0,jq1);
 177     qq02             = _mm_mul_ps(iq0,jq2);
 178     qq10             = _mm_mul_ps(iq1,jq0);
 179     qq11             = _mm_mul_ps(iq1,jq1);
 180     qq12             = _mm_mul_ps(iq1,jq2);
 181     qq20             = _mm_mul_ps(iq2,jq0);
 182     qq21             = _mm_mul_ps(iq2,jq1);
 183     qq22             = _mm_mul_ps(iq2,jq2);
 184
 185     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
 186     rcutoff_scalar   = fr->ic->rcoulomb;
 187     rcutoff          = _mm_set1_ps(rcutoff_scalar);
 188     rcutoff2         = _mm_mul_ps(rcutoff,rcutoff);
 189
 190     sh_vdw_invrcut6  = _mm_set1_ps(fr->ic->sh_invrc6);
 191     rvdw             = _mm_set1_ps(fr->ic->rvdw);
 192
 193     /* Avoid stupid compiler warnings */
 194     jnrA = jnrB = jnrC = jnrD = 0;
 195     j_coord_offsetA = 0;
 196     j_coord_offsetB = 0;
 197     j_coord_offsetC = 0;
 198     j_coord_offsetD = 0;
 199
 200     outeriter        = 0;
 201     inneriter        = 0;
 202
 203     for(iidx=0;iidx<4*DIM;iidx++)
 204     {
 205         scratch[iidx] = 0.0;
 206     }
 207
 208     /* Start outer loop over neighborlists */
 209     for(iidx=0; iidx<nri; iidx++)
 210     {
 211         /* Load shift vector for this list */
 212         i_shift_offset   = DIM*shiftidx[iidx];
 213
 214         /* Load limits for loop over neighbors */
 215         j_index_start    = jindex[iidx];
 216         j_index_end      = jindex[iidx+1];
 217
 218         /* Get outer coordinate index */
 219         inr              = iinr[iidx];
 220         i_coord_offset   = DIM*inr;
 221
 222         /* Load i particle coords and add shift vector */
 223         gmx_mm_load_shift_and_3rvec_broadcast_ps(shiftvec+i_shift_offset,x+i_coord_offset,
 224                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2);
 225
 226         fix0             = _mm_setzero_ps();
 227         fiy0             = _mm_setzero_ps();
 228         fiz0             = _mm_setzero_ps();
 229         fix1             = _mm_setzero_ps();
 230         fiy1             = _mm_setzero_ps();
 231         fiz1             = _mm_setzero_ps();
 232         fix2             = _mm_setzero_ps();
 233         fiy2             = _mm_setzero_ps();
 234         fiz2             = _mm_setzero_ps();
 235
 236         /* Reset potential sums */
 237         velecsum         = _mm_setzero_ps();
 238         vvdwsum          = _mm_setzero_ps();
 239
 240         /* Start inner kernel loop */
 241         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
 242         {
 243
 244             /* Get j neighbor index, and coordinate index */
 245             jnrA             = jjnr[jidx];
 246             jnrB             = jjnr[jidx+1];
 247             jnrC             = jjnr[jidx+2];
 248             jnrD             = jjnr[jidx+3];
 249             j_coord_offsetA  = DIM*jnrA;
 250             j_coord_offsetB  = DIM*jnrB;
 251             j_coord_offsetC  = DIM*jnrC;
 252             j_coord_offsetD  = DIM*jnrD;
 253
 254             /* load j atom coordinates */
 255             gmx_mm_load_3rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
 256                                               x+j_coord_offsetC,x+j_coord_offsetD,
 257                                               &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,&jy2,&jz2);
 258
 259             /* Calculate displacement vector */
 260             dx00             = _mm_sub_ps(ix0,jx0);
 261             dy00             = _mm_sub_ps(iy0,jy0);
 262             dz00             = _mm_sub_ps(iz0,jz0);
 263             dx01             = _mm_sub_ps(ix0,jx1);
 264             dy01             = _mm_sub_ps(iy0,jy1);
 265             dz01             = _mm_sub_ps(iz0,jz1);
 266             dx02             = _mm_sub_ps(ix0,jx2);
 267             dy02             = _mm_sub_ps(iy0,jy2);
 268             dz02             = _mm_sub_ps(iz0,jz2);
 269             dx10             = _mm_sub_ps(ix1,jx0);
 270             dy10             = _mm_sub_ps(iy1,jy0);
 271             dz10             = _mm_sub_ps(iz1,jz0);
 272             dx11             = _mm_sub_ps(ix1,jx1);
 273             dy11             = _mm_sub_ps(iy1,jy1);
 274             dz11             = _mm_sub_ps(iz1,jz1);
 275             dx12             = _mm_sub_ps(ix1,jx2);
 276             dy12             = _mm_sub_ps(iy1,jy2);
 277             dz12             = _mm_sub_ps(iz1,jz2);
 278             dx20             = _mm_sub_ps(ix2,jx0);
 279             dy20             = _mm_sub_ps(iy2,jy0);
 280             dz20             = _mm_sub_ps(iz2,jz0);
 281             dx21             = _mm_sub_ps(ix2,jx1);
 282             dy21             = _mm_sub_ps(iy2,jy1);
 283             dz21             = _mm_sub_ps(iz2,jz1);
 284             dx22             = _mm_sub_ps(ix2,jx2);
 285             dy22             = _mm_sub_ps(iy2,jy2);
 286             dz22             = _mm_sub_ps(iz2,jz2);
 287
 288             /* Calculate squared distance and things based on it */
 289             rsq00            = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
 290             rsq01            = gmx_mm_calc_rsq_ps(dx01,dy01,dz01);
 291             rsq02            = gmx_mm_calc_rsq_ps(dx02,dy02,dz02);
 292             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
 293             rsq11            = gmx_mm_calc_rsq_ps(dx11,dy11,dz11);
 294             rsq12            = gmx_mm_calc_rsq_ps(dx12,dy12,dz12);
 295             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
 296             rsq21            = gmx_mm_calc_rsq_ps(dx21,dy21,dz21);
 297             rsq22            = gmx_mm_calc_rsq_ps(dx22,dy22,dz22);
 298
 299             rinv00           = avx128fma_invsqrt_f(rsq00);
 300             rinv01           = avx128fma_invsqrt_f(rsq01);
 301             rinv02           = avx128fma_invsqrt_f(rsq02);
 302             rinv10           = avx128fma_invsqrt_f(rsq10);
 303             rinv11           = avx128fma_invsqrt_f(rsq11);
 304             rinv12           = avx128fma_invsqrt_f(rsq12);
 305             rinv20           = avx128fma_invsqrt_f(rsq20);
 306             rinv21           = avx128fma_invsqrt_f(rsq21);
 307             rinv22           = avx128fma_invsqrt_f(rsq22);
 308
 309             rinvsq00         = _mm_mul_ps(rinv00,rinv00);
 310             rinvsq01         = _mm_mul_ps(rinv01,rinv01);
 311             rinvsq02         = _mm_mul_ps(rinv02,rinv02);
 312             rinvsq10         = _mm_mul_ps(rinv10,rinv10);
 313             rinvsq11         = _mm_mul_ps(rinv11,rinv11);
 314             rinvsq12         = _mm_mul_ps(rinv12,rinv12);
 315             rinvsq20         = _mm_mul_ps(rinv20,rinv20);
 316             rinvsq21         = _mm_mul_ps(rinv21,rinv21);
 317             rinvsq22         = _mm_mul_ps(rinv22,rinv22);
 318
 319             fjx0             = _mm_setzero_ps();
 320             fjy0             = _mm_setzero_ps();
 321             fjz0             = _mm_setzero_ps();
 322             fjx1             = _mm_setzero_ps();
 323             fjy1             = _mm_setzero_ps();
 324             fjz1             = _mm_setzero_ps();
 325             fjx2             = _mm_setzero_ps();
 326             fjy2             = _mm_setzero_ps();
 327             fjz2             = _mm_setzero_ps();
 328
 329             /**************************
 330              * CALCULATE INTERACTIONS *
 331              **************************/
 332
 333             if (gmx_mm_any_lt(rsq00,rcutoff2))
 334             {
 335
 336             r00              = _mm_mul_ps(rsq00,rinv00);
 337
 338             /* EWALD ELECTROSTATICS */
 339
 340             /* Analytical PME correction */
 341             zeta2            = _mm_mul_ps(beta2,rsq00);
 342             rinv3            = _mm_mul_ps(rinvsq00,rinv00);
 343             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
 344             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
 345             felec            = _mm_mul_ps(qq00,felec);
 346             pmecorrV         = avx128fma_pmecorrV_f(zeta2);
 347             velec            = _mm_nmacc_ps(pmecorrV,beta,_mm_sub_ps(rinv00,sh_ewald));
 348             velec            = _mm_mul_ps(qq00,velec);
 349
 350             /* Analytical LJ-PME */
 351             rinvsix          = _mm_mul_ps(_mm_mul_ps(rinvsq00,rinvsq00),rinvsq00);
 352             ewcljrsq         = _mm_mul_ps(ewclj2,rsq00);
 353             ewclj6           = _mm_mul_ps(ewclj2,_mm_mul_ps(ewclj2,ewclj2));
 354             exponent         = avx128fma_exp_f(ewcljrsq);
 355             /* poly = exp(-(beta*r)^2) * (1 + (beta*r)^2 + (beta*r)^4 /2) */
 356             poly             = _mm_mul_ps(exponent,_mm_macc_ps(_mm_mul_ps(ewcljrsq,ewcljrsq),one_half,_mm_sub_ps(one,ewcljrsq)));
 357             /* vvdw6 = [C6 - C6grid * (1-poly)]/r6 */
 358             vvdw6            = _mm_mul_ps(_mm_macc_ps(-c6grid_00,_mm_sub_ps(one,poly),c6_00),rinvsix);
 359             vvdw12           = _mm_mul_ps(c12_00,_mm_mul_ps(rinvsix,rinvsix));
 360             vvdw             = _mm_msub_ps(_mm_nmacc_ps(c12_00,_mm_mul_ps(sh_vdw_invrcut6,sh_vdw_invrcut6),vvdw12),one_twelfth,
 361                                           _mm_mul_ps(_mm_sub_ps(vvdw6,_mm_macc_ps(c6grid_00,sh_lj_ewald,_mm_mul_ps(c6_00,sh_vdw_invrcut6))),one_sixth));
 362             /* fvdw = vvdw12/r - (vvdw6/r + (C6grid * exponent * beta^6)/r) */
 363             fvdw             = _mm_mul_ps(_mm_add_ps(vvdw12,_mm_msub_ps(_mm_mul_ps(c6grid_00,one_sixth),_mm_mul_ps(exponent,ewclj6),vvdw6)),rinvsq00);
 364
 365             cutoff_mask      = _mm_cmplt_ps(rsq00,rcutoff2);
 366
 367             /* Update potential sum for this i atom from the interaction with this j atom. */
 368             velec            = _mm_and_ps(velec,cutoff_mask);
 369             velecsum         = _mm_add_ps(velecsum,velec);
 370             vvdw             = _mm_and_ps(vvdw,cutoff_mask);
 371             vvdwsum          = _mm_add_ps(vvdwsum,vvdw);
 372
 373             fscal            = _mm_add_ps(felec,fvdw);
 374
 375             fscal            = _mm_and_ps(fscal,cutoff_mask);
 376
 377              /* Update vectorial force */
 378             fix0             = _mm_macc_ps(dx00,fscal,fix0);
 379             fiy0             = _mm_macc_ps(dy00,fscal,fiy0);
 380             fiz0             = _mm_macc_ps(dz00,fscal,fiz0);
 381
 382             fjx0             = _mm_macc_ps(dx00,fscal,fjx0);
 383             fjy0             = _mm_macc_ps(dy00,fscal,fjy0);
 384             fjz0             = _mm_macc_ps(dz00,fscal,fjz0);
 385
 386             }
 387
 388             /**************************
 389              * CALCULATE INTERACTIONS *
 390              **************************/
 391
 392             if (gmx_mm_any_lt(rsq01,rcutoff2))
 393             {
 394
 395             r01              = _mm_mul_ps(rsq01,rinv01);
 396
 397             /* EWALD ELECTROSTATICS */
 398
 399             /* Analytical PME correction */
 400             zeta2            = _mm_mul_ps(beta2,rsq01);
 401             rinv3            = _mm_mul_ps(rinvsq01,rinv01);
 402             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
 403             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
 404             felec            = _mm_mul_ps(qq01,felec);
 405             pmecorrV         = avx128fma_pmecorrV_f(zeta2);
 406             velec            = _mm_nmacc_ps(pmecorrV,beta,_mm_sub_ps(rinv01,sh_ewald));
 407             velec            = _mm_mul_ps(qq01,velec);
 408
 409             cutoff_mask      = _mm_cmplt_ps(rsq01,rcutoff2);
 410
 411             /* Update potential sum for this i atom from the interaction with this j atom. */
 412             velec            = _mm_and_ps(velec,cutoff_mask);
 413             velecsum         = _mm_add_ps(velecsum,velec);
 414
 415             fscal            = felec;
 416
 417             fscal            = _mm_and_ps(fscal,cutoff_mask);
 418
 419              /* Update vectorial force */
 420             fix0             = _mm_macc_ps(dx01,fscal,fix0);
 421             fiy0             = _mm_macc_ps(dy01,fscal,fiy0);
 422             fiz0             = _mm_macc_ps(dz01,fscal,fiz0);
 423
 424             fjx1             = _mm_macc_ps(dx01,fscal,fjx1);
 425             fjy1             = _mm_macc_ps(dy01,fscal,fjy1);
 426             fjz1             = _mm_macc_ps(dz01,fscal,fjz1);
 427
 428             }
 429
 430             /**************************
 431              * CALCULATE INTERACTIONS *
 432              **************************/
 433
 434             if (gmx_mm_any_lt(rsq02,rcutoff2))
 435             {
 436
 437             r02              = _mm_mul_ps(rsq02,rinv02);
 438
 439             /* EWALD ELECTROSTATICS */
 440
 441             /* Analytical PME correction */
 442             zeta2            = _mm_mul_ps(beta2,rsq02);
 443             rinv3            = _mm_mul_ps(rinvsq02,rinv02);
 444             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
 445             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
 446             felec            = _mm_mul_ps(qq02,felec);
 447             pmecorrV         = avx128fma_pmecorrV_f(zeta2);
 448             velec            = _mm_nmacc_ps(pmecorrV,beta,_mm_sub_ps(rinv02,sh_ewald));
 449             velec            = _mm_mul_ps(qq02,velec);
 450
 451             cutoff_mask      = _mm_cmplt_ps(rsq02,rcutoff2);
 452
 453             /* Update potential sum for this i atom from the interaction with this j atom. */
 454             velec            = _mm_and_ps(velec,cutoff_mask);
 455             velecsum         = _mm_add_ps(velecsum,velec);
 456
 457             fscal            = felec;
 458
 459             fscal            = _mm_and_ps(fscal,cutoff_mask);
 460
 461              /* Update vectorial force */
 462             fix0             = _mm_macc_ps(dx02,fscal,fix0);
 463             fiy0             = _mm_macc_ps(dy02,fscal,fiy0);
 464             fiz0             = _mm_macc_ps(dz02,fscal,fiz0);
 465
 466             fjx2             = _mm_macc_ps(dx02,fscal,fjx2);
 467             fjy2             = _mm_macc_ps(dy02,fscal,fjy2);
 468             fjz2             = _mm_macc_ps(dz02,fscal,fjz2);
 469
 470             }
 471
 472             /**************************
 473              * CALCULATE INTERACTIONS *
 474              **************************/
 475
 476             if (gmx_mm_any_lt(rsq10,rcutoff2))
 477             {
 478
 479             r10              = _mm_mul_ps(rsq10,rinv10);
 480
 481             /* EWALD ELECTROSTATICS */
 482
 483             /* Analytical PME correction */
 484             zeta2            = _mm_mul_ps(beta2,rsq10);
 485             rinv3            = _mm_mul_ps(rinvsq10,rinv10);
 486             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
 487             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
 488             felec            = _mm_mul_ps(qq10,felec);
 489             pmecorrV         = avx128fma_pmecorrV_f(zeta2);
 490             velec            = _mm_nmacc_ps(pmecorrV,beta,_mm_sub_ps(rinv10,sh_ewald));
 491             velec            = _mm_mul_ps(qq10,velec);
 492
 493             cutoff_mask      = _mm_cmplt_ps(rsq10,rcutoff2);
 494
 495             /* Update potential sum for this i atom from the interaction with this j atom. */
 496             velec            = _mm_and_ps(velec,cutoff_mask);
 497             velecsum         = _mm_add_ps(velecsum,velec);
 498
 499             fscal            = felec;
 500
 501             fscal            = _mm_and_ps(fscal,cutoff_mask);
 502
 503              /* Update vectorial force */
 504             fix1             = _mm_macc_ps(dx10,fscal,fix1);
 505             fiy1             = _mm_macc_ps(dy10,fscal,fiy1);
 506             fiz1             = _mm_macc_ps(dz10,fscal,fiz1);
 507
 508             fjx0             = _mm_macc_ps(dx10,fscal,fjx0);
 509             fjy0             = _mm_macc_ps(dy10,fscal,fjy0);
 510             fjz0             = _mm_macc_ps(dz10,fscal,fjz0);
 511
 512             }
 513
 514             /**************************
 515              * CALCULATE INTERACTIONS *
 516              **************************/
 517
 518             if (gmx_mm_any_lt(rsq11,rcutoff2))
 519             {
 520
 521             r11              = _mm_mul_ps(rsq11,rinv11);
 522
 523             /* EWALD ELECTROSTATICS */
 524
 525             /* Analytical PME correction */
 526             zeta2            = _mm_mul_ps(beta2,rsq11);
 527             rinv3            = _mm_mul_ps(rinvsq11,rinv11);
 528             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
 529             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
 530             felec            = _mm_mul_ps(qq11,felec);
 531             pmecorrV         = avx128fma_pmecorrV_f(zeta2);
 532             velec            = _mm_nmacc_ps(pmecorrV,beta,_mm_sub_ps(rinv11,sh_ewald));
 533             velec            = _mm_mul_ps(qq11,velec);
 534
 535             cutoff_mask      = _mm_cmplt_ps(rsq11,rcutoff2);
 536
 537             /* Update potential sum for this i atom from the interaction with this j atom. */
 538             velec            = _mm_and_ps(velec,cutoff_mask);
 539             velecsum         = _mm_add_ps(velecsum,velec);
 540
 541             fscal            = felec;
 542
 543             fscal            = _mm_and_ps(fscal,cutoff_mask);
 544
 545              /* Update vectorial force */
 546             fix1             = _mm_macc_ps(dx11,fscal,fix1);
 547             fiy1             = _mm_macc_ps(dy11,fscal,fiy1);
 548             fiz1             = _mm_macc_ps(dz11,fscal,fiz1);
 549
 550             fjx1             = _mm_macc_ps(dx11,fscal,fjx1);
 551             fjy1             = _mm_macc_ps(dy11,fscal,fjy1);
 552             fjz1             = _mm_macc_ps(dz11,fscal,fjz1);
 553
 554             }
 555
 556             /**************************
 557              * CALCULATE INTERACTIONS *
 558              **************************/
 559
 560             if (gmx_mm_any_lt(rsq12,rcutoff2))
 561             {
 562
 563             r12              = _mm_mul_ps(rsq12,rinv12);
 564
 565             /* EWALD ELECTROSTATICS */
 566
 567             /* Analytical PME correction */
 568             zeta2            = _mm_mul_ps(beta2,rsq12);
 569             rinv3            = _mm_mul_ps(rinvsq12,rinv12);
 570             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
 571             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
 572             felec            = _mm_mul_ps(qq12,felec);
 573             pmecorrV         = avx128fma_pmecorrV_f(zeta2);
 574             velec            = _mm_nmacc_ps(pmecorrV,beta,_mm_sub_ps(rinv12,sh_ewald));
 575             velec            = _mm_mul_ps(qq12,velec);
 576
 577             cutoff_mask      = _mm_cmplt_ps(rsq12,rcutoff2);
 578
 579             /* Update potential sum for this i atom from the interaction with this j atom. */
 580             velec            = _mm_and_ps(velec,cutoff_mask);
 581             velecsum         = _mm_add_ps(velecsum,velec);
 582
 583             fscal            = felec;
 584
 585             fscal            = _mm_and_ps(fscal,cutoff_mask);
 586
 587              /* Update vectorial force */
 588             fix1             = _mm_macc_ps(dx12,fscal,fix1);
 589             fiy1             = _mm_macc_ps(dy12,fscal,fiy1);
 590             fiz1             = _mm_macc_ps(dz12,fscal,fiz1);
 591
 592             fjx2             = _mm_macc_ps(dx12,fscal,fjx2);
 593             fjy2             = _mm_macc_ps(dy12,fscal,fjy2);
 594             fjz2             = _mm_macc_ps(dz12,fscal,fjz2);
 595
 596             }
 597
 598             /**************************
 599              * CALCULATE INTERACTIONS *
 600              **************************/
 601
 602             if (gmx_mm_any_lt(rsq20,rcutoff2))
 603             {
 604
 605             r20              = _mm_mul_ps(rsq20,rinv20);
 606
 607             /* EWALD ELECTROSTATICS */
 608
 609             /* Analytical PME correction */
 610             zeta2            = _mm_mul_ps(beta2,rsq20);
 611             rinv3            = _mm_mul_ps(rinvsq20,rinv20);
 612             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
 613             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
 614             felec            = _mm_mul_ps(qq20,felec);
 615             pmecorrV         = avx128fma_pmecorrV_f(zeta2);
 616             velec            = _mm_nmacc_ps(pmecorrV,beta,_mm_sub_ps(rinv20,sh_ewald));
 617             velec            = _mm_mul_ps(qq20,velec);
 618
 619             cutoff_mask      = _mm_cmplt_ps(rsq20,rcutoff2);
 620
 621             /* Update potential sum for this i atom from the interaction with this j atom. */
 622             velec            = _mm_and_ps(velec,cutoff_mask);
 623             velecsum         = _mm_add_ps(velecsum,velec);
 624
 625             fscal            = felec;
 626
 627             fscal            = _mm_and_ps(fscal,cutoff_mask);
 628
 629              /* Update vectorial force */
 630             fix2             = _mm_macc_ps(dx20,fscal,fix2);
 631             fiy2             = _mm_macc_ps(dy20,fscal,fiy2);
 632             fiz2             = _mm_macc_ps(dz20,fscal,fiz2);
 633
 634             fjx0             = _mm_macc_ps(dx20,fscal,fjx0);
 635             fjy0             = _mm_macc_ps(dy20,fscal,fjy0);
 636             fjz0             = _mm_macc_ps(dz20,fscal,fjz0);
 637
 638             }
 639
 640             /**************************
 641              * CALCULATE INTERACTIONS *
 642              **************************/
 643
 644             if (gmx_mm_any_lt(rsq21,rcutoff2))
 645             {
 646
 647             r21              = _mm_mul_ps(rsq21,rinv21);
 648
 649             /* EWALD ELECTROSTATICS */
 650
 651             /* Analytical PME correction */
 652             zeta2            = _mm_mul_ps(beta2,rsq21);
 653             rinv3            = _mm_mul_ps(rinvsq21,rinv21);
 654             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
 655             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
 656             felec            = _mm_mul_ps(qq21,felec);
 657             pmecorrV         = avx128fma_pmecorrV_f(zeta2);
 658             velec            = _mm_nmacc_ps(pmecorrV,beta,_mm_sub_ps(rinv21,sh_ewald));
 659             velec            = _mm_mul_ps(qq21,velec);
 660
 661             cutoff_mask      = _mm_cmplt_ps(rsq21,rcutoff2);
 662
 663             /* Update potential sum for this i atom from the interaction with this j atom. */
 664             velec            = _mm_and_ps(velec,cutoff_mask);
 665             velecsum         = _mm_add_ps(velecsum,velec);
 666
 667             fscal            = felec;
 668
 669             fscal            = _mm_and_ps(fscal,cutoff_mask);
 670
 671              /* Update vectorial force */
 672             fix2             = _mm_macc_ps(dx21,fscal,fix2);
 673             fiy2             = _mm_macc_ps(dy21,fscal,fiy2);
 674             fiz2             = _mm_macc_ps(dz21,fscal,fiz2);
 675
 676             fjx1             = _mm_macc_ps(dx21,fscal,fjx1);
 677             fjy1             = _mm_macc_ps(dy21,fscal,fjy1);
 678             fjz1             = _mm_macc_ps(dz21,fscal,fjz1);
 679
 680             }
 681
 682             /**************************
 683              * CALCULATE INTERACTIONS *
 684              **************************/
 685
 686             if (gmx_mm_any_lt(rsq22,rcutoff2))
 687             {
 688
 689             r22              = _mm_mul_ps(rsq22,rinv22);
 690
 691             /* EWALD ELECTROSTATICS */
 692
 693             /* Analytical PME correction */
 694             zeta2            = _mm_mul_ps(beta2,rsq22);
 695             rinv3            = _mm_mul_ps(rinvsq22,rinv22);
 696             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
 697             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
 698             felec            = _mm_mul_ps(qq22,felec);
 699             pmecorrV         = avx128fma_pmecorrV_f(zeta2);
 700             velec            = _mm_nmacc_ps(pmecorrV,beta,_mm_sub_ps(rinv22,sh_ewald));
 701             velec            = _mm_mul_ps(qq22,velec);
 702
 703             cutoff_mask      = _mm_cmplt_ps(rsq22,rcutoff2);
 704
 705             /* Update potential sum for this i atom from the interaction with this j atom. */
 706             velec            = _mm_and_ps(velec,cutoff_mask);
 707             velecsum         = _mm_add_ps(velecsum,velec);
 708
 709             fscal            = felec;
 710
 711             fscal            = _mm_and_ps(fscal,cutoff_mask);
 712
 713              /* Update vectorial force */
 714             fix2             = _mm_macc_ps(dx22,fscal,fix2);
 715             fiy2             = _mm_macc_ps(dy22,fscal,fiy2);
 716             fiz2             = _mm_macc_ps(dz22,fscal,fiz2);
 717
 718             fjx2             = _mm_macc_ps(dx22,fscal,fjx2);
 719             fjy2             = _mm_macc_ps(dy22,fscal,fjy2);
 720             fjz2             = _mm_macc_ps(dz22,fscal,fjz2);
 721
 722             }
 723
 724             fjptrA             = f+j_coord_offsetA;
 725             fjptrB             = f+j_coord_offsetB;
 726             fjptrC             = f+j_coord_offsetC;
 727             fjptrD             = f+j_coord_offsetD;
 728
 729             gmx_mm_decrement_3rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,
 730                                                    fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,fjx2,fjy2,fjz2);
 731
 732             /* Inner loop uses 327 flops */
 733         }
 734
 735         if(jidx<j_index_end)
 736         {
 737
 738             /* Get j neighbor index, and coordinate index */
 739             jnrlistA         = jjnr[jidx];
 740             jnrlistB         = jjnr[jidx+1];
 741             jnrlistC         = jjnr[jidx+2];
 742             jnrlistD         = jjnr[jidx+3];
 743             /* Sign of each element will be negative for non-real atoms.
 744              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
 745              * so use it as val = _mm_andnot_ps(mask,val) to clear dummy entries.
 746              */
 747             dummy_mask = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
 748             jnrA       = (jnrlistA>=0) ? jnrlistA : 0;
 749             jnrB       = (jnrlistB>=0) ? jnrlistB : 0;
 750             jnrC       = (jnrlistC>=0) ? jnrlistC : 0;
 751             jnrD       = (jnrlistD>=0) ? jnrlistD : 0;
 752             j_coord_offsetA  = DIM*jnrA;
 753             j_coord_offsetB  = DIM*jnrB;
 754             j_coord_offsetC  = DIM*jnrC;
 755             j_coord_offsetD  = DIM*jnrD;
 756
 757             /* load j atom coordinates */
 758             gmx_mm_load_3rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
 759                                               x+j_coord_offsetC,x+j_coord_offsetD,
 760                                               &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,&jy2,&jz2);
 761
 762             /* Calculate displacement vector */
 763             dx00             = _mm_sub_ps(ix0,jx0);
 764             dy00             = _mm_sub_ps(iy0,jy0);
 765             dz00             = _mm_sub_ps(iz0,jz0);
 766             dx01             = _mm_sub_ps(ix0,jx1);
 767             dy01             = _mm_sub_ps(iy0,jy1);
 768             dz01             = _mm_sub_ps(iz0,jz1);
 769             dx02             = _mm_sub_ps(ix0,jx2);
 770             dy02             = _mm_sub_ps(iy0,jy2);
 771             dz02             = _mm_sub_ps(iz0,jz2);
 772             dx10             = _mm_sub_ps(ix1,jx0);
 773             dy10             = _mm_sub_ps(iy1,jy0);
 774             dz10             = _mm_sub_ps(iz1,jz0);
 775             dx11             = _mm_sub_ps(ix1,jx1);
 776             dy11             = _mm_sub_ps(iy1,jy1);
 777             dz11             = _mm_sub_ps(iz1,jz1);
 778             dx12             = _mm_sub_ps(ix1,jx2);
 779             dy12             = _mm_sub_ps(iy1,jy2);
 780             dz12             = _mm_sub_ps(iz1,jz2);
 781             dx20             = _mm_sub_ps(ix2,jx0);
 782             dy20             = _mm_sub_ps(iy2,jy0);
 783             dz20             = _mm_sub_ps(iz2,jz0);
 784             dx21             = _mm_sub_ps(ix2,jx1);
 785             dy21             = _mm_sub_ps(iy2,jy1);
 786             dz21             = _mm_sub_ps(iz2,jz1);
 787             dx22             = _mm_sub_ps(ix2,jx2);
 788             dy22             = _mm_sub_ps(iy2,jy2);
 789             dz22             = _mm_sub_ps(iz2,jz2);
 790
 791             /* Calculate squared distance and things based on it */
 792             rsq00            = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
 793             rsq01            = gmx_mm_calc_rsq_ps(dx01,dy01,dz01);
 794             rsq02            = gmx_mm_calc_rsq_ps(dx02,dy02,dz02);
 795             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
 796             rsq11            = gmx_mm_calc_rsq_ps(dx11,dy11,dz11);
 797             rsq12            = gmx_mm_calc_rsq_ps(dx12,dy12,dz12);
 798             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
 799             rsq21            = gmx_mm_calc_rsq_ps(dx21,dy21,dz21);
 800             rsq22            = gmx_mm_calc_rsq_ps(dx22,dy22,dz22);
 801
 802             rinv00           = avx128fma_invsqrt_f(rsq00);
 803             rinv01           = avx128fma_invsqrt_f(rsq01);
 804             rinv02           = avx128fma_invsqrt_f(rsq02);
 805             rinv10           = avx128fma_invsqrt_f(rsq10);
 806             rinv11           = avx128fma_invsqrt_f(rsq11);
 807             rinv12           = avx128fma_invsqrt_f(rsq12);
 808             rinv20           = avx128fma_invsqrt_f(rsq20);
 809             rinv21           = avx128fma_invsqrt_f(rsq21);
 810             rinv22           = avx128fma_invsqrt_f(rsq22);
 811
 812             rinvsq00         = _mm_mul_ps(rinv00,rinv00);
 813             rinvsq01         = _mm_mul_ps(rinv01,rinv01);
 814             rinvsq02         = _mm_mul_ps(rinv02,rinv02);
 815             rinvsq10         = _mm_mul_ps(rinv10,rinv10);
 816             rinvsq11         = _mm_mul_ps(rinv11,rinv11);
 817             rinvsq12         = _mm_mul_ps(rinv12,rinv12);
 818             rinvsq20         = _mm_mul_ps(rinv20,rinv20);
 819             rinvsq21         = _mm_mul_ps(rinv21,rinv21);
 820             rinvsq22         = _mm_mul_ps(rinv22,rinv22);
 821
 822             fjx0             = _mm_setzero_ps();
 823             fjy0             = _mm_setzero_ps();
 824             fjz0             = _mm_setzero_ps();
 825             fjx1             = _mm_setzero_ps();
 826             fjy1             = _mm_setzero_ps();
 827             fjz1             = _mm_setzero_ps();
 828             fjx2             = _mm_setzero_ps();
 829             fjy2             = _mm_setzero_ps();
 830             fjz2             = _mm_setzero_ps();
 831
 832             /**************************
 833              * CALCULATE INTERACTIONS *
 834              **************************/
 835
 836             if (gmx_mm_any_lt(rsq00,rcutoff2))
 837             {
 838
 839             r00              = _mm_mul_ps(rsq00,rinv00);
 840             r00              = _mm_andnot_ps(dummy_mask,r00);
 841
 842             /* EWALD ELECTROSTATICS */
 843
 844             /* Analytical PME correction */
 845             zeta2            = _mm_mul_ps(beta2,rsq00);
 846             rinv3            = _mm_mul_ps(rinvsq00,rinv00);
 847             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
 848             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
 849             felec            = _mm_mul_ps(qq00,felec);
 850             pmecorrV         = avx128fma_pmecorrV_f(zeta2);
 851             velec            = _mm_nmacc_ps(pmecorrV,beta,_mm_sub_ps(rinv00,sh_ewald));
 852             velec            = _mm_mul_ps(qq00,velec);
 853
 854             /* Analytical LJ-PME */
 855             rinvsix          = _mm_mul_ps(_mm_mul_ps(rinvsq00,rinvsq00),rinvsq00);
 856             ewcljrsq         = _mm_mul_ps(ewclj2,rsq00);
 857             ewclj6           = _mm_mul_ps(ewclj2,_mm_mul_ps(ewclj2,ewclj2));
 858             exponent         = avx128fma_exp_f(ewcljrsq);
 859             /* poly = exp(-(beta*r)^2) * (1 + (beta*r)^2 + (beta*r)^4 /2) */
 860             poly             = _mm_mul_ps(exponent,_mm_macc_ps(_mm_mul_ps(ewcljrsq,ewcljrsq),one_half,_mm_sub_ps(one,ewcljrsq)));
 861             /* vvdw6 = [C6 - C6grid * (1-poly)]/r6 */
 862             vvdw6            = _mm_mul_ps(_mm_macc_ps(-c6grid_00,_mm_sub_ps(one,poly),c6_00),rinvsix);
 863             vvdw12           = _mm_mul_ps(c12_00,_mm_mul_ps(rinvsix,rinvsix));
 864             vvdw             = _mm_msub_ps(_mm_nmacc_ps(c12_00,_mm_mul_ps(sh_vdw_invrcut6,sh_vdw_invrcut6),vvdw12),one_twelfth,
 865                                           _mm_mul_ps(_mm_sub_ps(vvdw6,_mm_macc_ps(c6grid_00,sh_lj_ewald,_mm_mul_ps(c6_00,sh_vdw_invrcut6))),one_sixth));
 866             /* fvdw = vvdw12/r - (vvdw6/r + (C6grid * exponent * beta^6)/r) */
 867             fvdw             = _mm_mul_ps(_mm_add_ps(vvdw12,_mm_msub_ps(_mm_mul_ps(c6grid_00,one_sixth),_mm_mul_ps(exponent,ewclj6),vvdw6)),rinvsq00);
 868
 869             cutoff_mask      = _mm_cmplt_ps(rsq00,rcutoff2);
 870
 871             /* Update potential sum for this i atom from the interaction with this j atom. */
 872             velec            = _mm_and_ps(velec,cutoff_mask);
 873             velec            = _mm_andnot_ps(dummy_mask,velec);
 874             velecsum         = _mm_add_ps(velecsum,velec);
 875             vvdw             = _mm_and_ps(vvdw,cutoff_mask);
 876             vvdw             = _mm_andnot_ps(dummy_mask,vvdw);
 877             vvdwsum          = _mm_add_ps(vvdwsum,vvdw);
 878
 879             fscal            = _mm_add_ps(felec,fvdw);
 880
 881             fscal            = _mm_and_ps(fscal,cutoff_mask);
 882
 883             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 884
 885              /* Update vectorial force */
 886             fix0             = _mm_macc_ps(dx00,fscal,fix0);
 887             fiy0             = _mm_macc_ps(dy00,fscal,fiy0);
 888             fiz0             = _mm_macc_ps(dz00,fscal,fiz0);
 889
 890             fjx0             = _mm_macc_ps(dx00,fscal,fjx0);
 891             fjy0             = _mm_macc_ps(dy00,fscal,fjy0);
 892             fjz0             = _mm_macc_ps(dz00,fscal,fjz0);
 893
 894             }
 895
 896             /**************************
 897              * CALCULATE INTERACTIONS *
 898              **************************/
 899
 900             if (gmx_mm_any_lt(rsq01,rcutoff2))
 901             {
 902
 903             r01              = _mm_mul_ps(rsq01,rinv01);
 904             r01              = _mm_andnot_ps(dummy_mask,r01);
 905
 906             /* EWALD ELECTROSTATICS */
 907
 908             /* Analytical PME correction */
 909             zeta2            = _mm_mul_ps(beta2,rsq01);
 910             rinv3            = _mm_mul_ps(rinvsq01,rinv01);
 911             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
 912             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
 913             felec            = _mm_mul_ps(qq01,felec);
 914             pmecorrV         = avx128fma_pmecorrV_f(zeta2);
 915             velec            = _mm_nmacc_ps(pmecorrV,beta,_mm_sub_ps(rinv01,sh_ewald));
 916             velec            = _mm_mul_ps(qq01,velec);
 917
 918             cutoff_mask      = _mm_cmplt_ps(rsq01,rcutoff2);
 919
 920             /* Update potential sum for this i atom from the interaction with this j atom. */
 921             velec            = _mm_and_ps(velec,cutoff_mask);
 922             velec            = _mm_andnot_ps(dummy_mask,velec);
 923             velecsum         = _mm_add_ps(velecsum,velec);
 924
 925             fscal            = felec;
 926
 927             fscal            = _mm_and_ps(fscal,cutoff_mask);
 928
 929             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 930
 931              /* Update vectorial force */
 932             fix0             = _mm_macc_ps(dx01,fscal,fix0);
 933             fiy0             = _mm_macc_ps(dy01,fscal,fiy0);
 934             fiz0             = _mm_macc_ps(dz01,fscal,fiz0);
 935
 936             fjx1             = _mm_macc_ps(dx01,fscal,fjx1);
 937             fjy1             = _mm_macc_ps(dy01,fscal,fjy1);
 938             fjz1             = _mm_macc_ps(dz01,fscal,fjz1);
 939
 940             }
 941
 942             /**************************
 943              * CALCULATE INTERACTIONS *
 944              **************************/
 945
 946             if (gmx_mm_any_lt(rsq02,rcutoff2))
 947             {
 948
 949             r02              = _mm_mul_ps(rsq02,rinv02);
 950             r02              = _mm_andnot_ps(dummy_mask,r02);
 951
 952             /* EWALD ELECTROSTATICS */
 953
 954             /* Analytical PME correction */
 955             zeta2            = _mm_mul_ps(beta2,rsq02);
 956             rinv3            = _mm_mul_ps(rinvsq02,rinv02);
 957             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
 958             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
 959             felec            = _mm_mul_ps(qq02,felec);
 960             pmecorrV         = avx128fma_pmecorrV_f(zeta2);
 961             velec            = _mm_nmacc_ps(pmecorrV,beta,_mm_sub_ps(rinv02,sh_ewald));
 962             velec            = _mm_mul_ps(qq02,velec);
 963
 964             cutoff_mask      = _mm_cmplt_ps(rsq02,rcutoff2);
 965
 966             /* Update potential sum for this i atom from the interaction with this j atom. */
 967             velec            = _mm_and_ps(velec,cutoff_mask);
 968             velec            = _mm_andnot_ps(dummy_mask,velec);
 969             velecsum         = _mm_add_ps(velecsum,velec);
 970
 971             fscal            = felec;
 972
 973             fscal            = _mm_and_ps(fscal,cutoff_mask);
 974
 975             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 976
 977              /* Update vectorial force */
 978             fix0             = _mm_macc_ps(dx02,fscal,fix0);
 979             fiy0             = _mm_macc_ps(dy02,fscal,fiy0);
 980             fiz0             = _mm_macc_ps(dz02,fscal,fiz0);
 981
 982             fjx2             = _mm_macc_ps(dx02,fscal,fjx2);
 983             fjy2             = _mm_macc_ps(dy02,fscal,fjy2);
 984             fjz2             = _mm_macc_ps(dz02,fscal,fjz2);
 985
 986             }
 987
 988             /**************************
 989              * CALCULATE INTERACTIONS *
 990              **************************/
 991
 992             if (gmx_mm_any_lt(rsq10,rcutoff2))
 993             {
 994
 995             r10              = _mm_mul_ps(rsq10,rinv10);
 996             r10              = _mm_andnot_ps(dummy_mask,r10);
 997
 998             /* EWALD ELECTROSTATICS */
 999
1000             /* Analytical PME correction */
1001             zeta2            = _mm_mul_ps(beta2,rsq10);
1002             rinv3            = _mm_mul_ps(rinvsq10,rinv10);
1003             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
1004             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
1005             felec            = _mm_mul_ps(qq10,felec);
1006             pmecorrV         = avx128fma_pmecorrV_f(zeta2);
1007             velec            = _mm_nmacc_ps(pmecorrV,beta,_mm_sub_ps(rinv10,sh_ewald));
1008             velec            = _mm_mul_ps(qq10,velec);
1009
1010             cutoff_mask      = _mm_cmplt_ps(rsq10,rcutoff2);
1011
1012             /* Update potential sum for this i atom from the interaction with this j atom. */
1013             velec            = _mm_and_ps(velec,cutoff_mask);
1014             velec            = _mm_andnot_ps(dummy_mask,velec);
1015             velecsum         = _mm_add_ps(velecsum,velec);
1016
1017             fscal            = felec;
1018
1019             fscal            = _mm_and_ps(fscal,cutoff_mask);
1020
1021             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1022
1023              /* Update vectorial force */
1024             fix1             = _mm_macc_ps(dx10,fscal,fix1);
1025             fiy1             = _mm_macc_ps(dy10,fscal,fiy1);
1026             fiz1             = _mm_macc_ps(dz10,fscal,fiz1);
1027
1028             fjx0             = _mm_macc_ps(dx10,fscal,fjx0);
1029             fjy0             = _mm_macc_ps(dy10,fscal,fjy0);
1030             fjz0             = _mm_macc_ps(dz10,fscal,fjz0);
1031
1032             }
1033
1034             /**************************
1035              * CALCULATE INTERACTIONS *
1036              **************************/
1037
1038             if (gmx_mm_any_lt(rsq11,rcutoff2))
1039             {
1040
1041             r11              = _mm_mul_ps(rsq11,rinv11);
1042             r11              = _mm_andnot_ps(dummy_mask,r11);
1043
1044             /* EWALD ELECTROSTATICS */
1045
1046             /* Analytical PME correction */
1047             zeta2            = _mm_mul_ps(beta2,rsq11);
1048             rinv3            = _mm_mul_ps(rinvsq11,rinv11);
1049             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
1050             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
1051             felec            = _mm_mul_ps(qq11,felec);
1052             pmecorrV         = avx128fma_pmecorrV_f(zeta2);
1053             velec            = _mm_nmacc_ps(pmecorrV,beta,_mm_sub_ps(rinv11,sh_ewald));
1054             velec            = _mm_mul_ps(qq11,velec);
1055
1056             cutoff_mask      = _mm_cmplt_ps(rsq11,rcutoff2);
1057
1058             /* Update potential sum for this i atom from the interaction with this j atom. */
1059             velec            = _mm_and_ps(velec,cutoff_mask);
1060             velec            = _mm_andnot_ps(dummy_mask,velec);
1061             velecsum         = _mm_add_ps(velecsum,velec);
1062
1063             fscal            = felec;
1064
1065             fscal            = _mm_and_ps(fscal,cutoff_mask);
1066
1067             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1068
1069              /* Update vectorial force */
1070             fix1             = _mm_macc_ps(dx11,fscal,fix1);
1071             fiy1             = _mm_macc_ps(dy11,fscal,fiy1);
1072             fiz1             = _mm_macc_ps(dz11,fscal,fiz1);
1073
1074             fjx1             = _mm_macc_ps(dx11,fscal,fjx1);
1075             fjy1             = _mm_macc_ps(dy11,fscal,fjy1);
1076             fjz1             = _mm_macc_ps(dz11,fscal,fjz1);
1077
1078             }
1079
1080             /**************************
1081              * CALCULATE INTERACTIONS *
1082              **************************/
1083
1084             if (gmx_mm_any_lt(rsq12,rcutoff2))
1085             {
1086
1087             r12              = _mm_mul_ps(rsq12,rinv12);
1088             r12              = _mm_andnot_ps(dummy_mask,r12);
1089
1090             /* EWALD ELECTROSTATICS */
1091
1092             /* Analytical PME correction */
1093             zeta2            = _mm_mul_ps(beta2,rsq12);
1094             rinv3            = _mm_mul_ps(rinvsq12,rinv12);
1095             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
1096             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
1097             felec            = _mm_mul_ps(qq12,felec);
1098             pmecorrV         = avx128fma_pmecorrV_f(zeta2);
1099             velec            = _mm_nmacc_ps(pmecorrV,beta,_mm_sub_ps(rinv12,sh_ewald));
1100             velec            = _mm_mul_ps(qq12,velec);
1101
1102             cutoff_mask      = _mm_cmplt_ps(rsq12,rcutoff2);
1103
1104             /* Update potential sum for this i atom from the interaction with this j atom. */
1105             velec            = _mm_and_ps(velec,cutoff_mask);
1106             velec            = _mm_andnot_ps(dummy_mask,velec);
1107             velecsum         = _mm_add_ps(velecsum,velec);
1108
1109             fscal            = felec;
1110
1111             fscal            = _mm_and_ps(fscal,cutoff_mask);
1112
1113             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1114
1115              /* Update vectorial force */
1116             fix1             = _mm_macc_ps(dx12,fscal,fix1);
1117             fiy1             = _mm_macc_ps(dy12,fscal,fiy1);
1118             fiz1             = _mm_macc_ps(dz12,fscal,fiz1);
1119
1120             fjx2             = _mm_macc_ps(dx12,fscal,fjx2);
1121             fjy2             = _mm_macc_ps(dy12,fscal,fjy2);
1122             fjz2             = _mm_macc_ps(dz12,fscal,fjz2);
1123
1124             }
1125
1126             /**************************
1127              * CALCULATE INTERACTIONS *
1128              **************************/
1129
1130             if (gmx_mm_any_lt(rsq20,rcutoff2))
1131             {
1132
1133             r20              = _mm_mul_ps(rsq20,rinv20);
1134             r20              = _mm_andnot_ps(dummy_mask,r20);
1135
1136             /* EWALD ELECTROSTATICS */
1137
1138             /* Analytical PME correction */
1139             zeta2            = _mm_mul_ps(beta2,rsq20);
1140             rinv3            = _mm_mul_ps(rinvsq20,rinv20);
1141             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
1142             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
1143             felec            = _mm_mul_ps(qq20,felec);
1144             pmecorrV         = avx128fma_pmecorrV_f(zeta2);
1145             velec            = _mm_nmacc_ps(pmecorrV,beta,_mm_sub_ps(rinv20,sh_ewald));
1146             velec            = _mm_mul_ps(qq20,velec);
1147
1148             cutoff_mask      = _mm_cmplt_ps(rsq20,rcutoff2);
1149
1150             /* Update potential sum for this i atom from the interaction with this j atom. */
1151             velec            = _mm_and_ps(velec,cutoff_mask);
1152             velec            = _mm_andnot_ps(dummy_mask,velec);
1153             velecsum         = _mm_add_ps(velecsum,velec);
1154
1155             fscal            = felec;
1156
1157             fscal            = _mm_and_ps(fscal,cutoff_mask);
1158
1159             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1160
1161              /* Update vectorial force */
1162             fix2             = _mm_macc_ps(dx20,fscal,fix2);
1163             fiy2             = _mm_macc_ps(dy20,fscal,fiy2);
1164             fiz2             = _mm_macc_ps(dz20,fscal,fiz2);
1165
1166             fjx0             = _mm_macc_ps(dx20,fscal,fjx0);
1167             fjy0             = _mm_macc_ps(dy20,fscal,fjy0);
1168             fjz0             = _mm_macc_ps(dz20,fscal,fjz0);
1169
1170             }
1171
1172             /**************************
1173              * CALCULATE INTERACTIONS *
1174              **************************/
1175
1176             if (gmx_mm_any_lt(rsq21,rcutoff2))
1177             {
1178
1179             r21              = _mm_mul_ps(rsq21,rinv21);
1180             r21              = _mm_andnot_ps(dummy_mask,r21);
1181
1182             /* EWALD ELECTROSTATICS */
1183
1184             /* Analytical PME correction */
1185             zeta2            = _mm_mul_ps(beta2,rsq21);
1186             rinv3            = _mm_mul_ps(rinvsq21,rinv21);
1187             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
1188             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
1189             felec            = _mm_mul_ps(qq21,felec);
1190             pmecorrV         = avx128fma_pmecorrV_f(zeta2);
1191             velec            = _mm_nmacc_ps(pmecorrV,beta,_mm_sub_ps(rinv21,sh_ewald));
1192             velec            = _mm_mul_ps(qq21,velec);
1193
1194             cutoff_mask      = _mm_cmplt_ps(rsq21,rcutoff2);
1195
1196             /* Update potential sum for this i atom from the interaction with this j atom. */
1197             velec            = _mm_and_ps(velec,cutoff_mask);
1198             velec            = _mm_andnot_ps(dummy_mask,velec);
1199             velecsum         = _mm_add_ps(velecsum,velec);
1200
1201             fscal            = felec;
1202
1203             fscal            = _mm_and_ps(fscal,cutoff_mask);
1204
1205             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1206
1207              /* Update vectorial force */
1208             fix2             = _mm_macc_ps(dx21,fscal,fix2);
1209             fiy2             = _mm_macc_ps(dy21,fscal,fiy2);
1210             fiz2             = _mm_macc_ps(dz21,fscal,fiz2);
1211
1212             fjx1             = _mm_macc_ps(dx21,fscal,fjx1);
1213             fjy1             = _mm_macc_ps(dy21,fscal,fjy1);
1214             fjz1             = _mm_macc_ps(dz21,fscal,fjz1);
1215
1216             }
1217
1218             /**************************
1219              * CALCULATE INTERACTIONS *
1220              **************************/
1221
1222             if (gmx_mm_any_lt(rsq22,rcutoff2))
1223             {
1224
1225             r22              = _mm_mul_ps(rsq22,rinv22);
1226             r22              = _mm_andnot_ps(dummy_mask,r22);
1227
1228             /* EWALD ELECTROSTATICS */
1229
1230             /* Analytical PME correction */
1231             zeta2            = _mm_mul_ps(beta2,rsq22);
1232             rinv3            = _mm_mul_ps(rinvsq22,rinv22);
1233             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
1234             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
1235             felec            = _mm_mul_ps(qq22,felec);
1236             pmecorrV         = avx128fma_pmecorrV_f(zeta2);
1237             velec            = _mm_nmacc_ps(pmecorrV,beta,_mm_sub_ps(rinv22,sh_ewald));
1238             velec            = _mm_mul_ps(qq22,velec);
1239
1240             cutoff_mask      = _mm_cmplt_ps(rsq22,rcutoff2);
1241
1242             /* Update potential sum for this i atom from the interaction with this j atom. */
1243             velec            = _mm_and_ps(velec,cutoff_mask);
1244             velec            = _mm_andnot_ps(dummy_mask,velec);
1245             velecsum         = _mm_add_ps(velecsum,velec);
1246
1247             fscal            = felec;
1248
1249             fscal            = _mm_and_ps(fscal,cutoff_mask);
1250
1251             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1252
1253              /* Update vectorial force */
1254             fix2             = _mm_macc_ps(dx22,fscal,fix2);
1255             fiy2             = _mm_macc_ps(dy22,fscal,fiy2);
1256             fiz2             = _mm_macc_ps(dz22,fscal,fiz2);
1257
1258             fjx2             = _mm_macc_ps(dx22,fscal,fjx2);
1259             fjy2             = _mm_macc_ps(dy22,fscal,fjy2);
1260             fjz2             = _mm_macc_ps(dz22,fscal,fjz2);
1261
1262             }
1263
1264             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
1265             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
1266             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
1267             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
1268
1269             gmx_mm_decrement_3rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,
1270                                                    fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,fjx2,fjy2,fjz2);
1271
1272             /* Inner loop uses 336 flops */
1273         }
1274
1275         /* End of innermost loop */
1276
1277         gmx_mm_update_iforce_3atom_swizzle_ps(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,
1278                                               f+i_coord_offset,fshift+i_shift_offset);
1279
1280         ggid                        = gid[iidx];
1281         /* Update potential energies */
1282         gmx_mm_update_1pot_ps(velecsum,kernel_data->energygrp_elec+ggid);
1283         gmx_mm_update_1pot_ps(vvdwsum,kernel_data->energygrp_vdw+ggid);
1284
1285         /* Increment number of inner iterations */
1286         inneriter                  += j_index_end - j_index_start;
1287
1288         /* Outer loop uses 20 flops */
1289     }
1290
1291     /* Increment number of outer iterations */
1292     outeriter        += nri;
1293
1294     /* Update outer/inner flops */
1295
1296     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W3W3_VF,outeriter*20 + inneriter*336);
1297 }
1298 /*
1299  * Gromacs nonbonded kernel:   nb_kernel_ElecEwSh_VdwLJEwSh_GeomW3W3_F_avx_128_fma_single
1300  * Electrostatics interaction: Ewald
1301  * VdW interaction:            LJEwald
1302  * Geometry:                   Water3-Water3
1303  * Calculate force/pot:        Force
1304  */
1305 void
1306 nb_kernel_ElecEwSh_VdwLJEwSh_GeomW3W3_F_avx_128_fma_single
1307                     (t_nblist                    * gmx_restrict       nlist,
1308                      rvec                        * gmx_restrict          xx,
1309                      rvec                        * gmx_restrict          ff,
1310                      struct t_forcerec           * gmx_restrict          fr,
1311                      t_mdatoms                   * gmx_restrict     mdatoms,
1312                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
1313                      t_nrnb                      * gmx_restrict        nrnb)
1314 {
1315     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
1316      * just 0 for non-waters.
1317      * Suffixes A,B,C,D refer to j loop unrolling done with AVX_128, e.g. for the four different
1318      * jnr indices corresponding to data put in the four positions in the SIMD register.
1319      */
1320     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
1321     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
1322     int              jnrA,jnrB,jnrC,jnrD;
1323     int              jnrlistA,jnrlistB,jnrlistC,jnrlistD;
1324     int              j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
1325     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
1326     real             rcutoff_scalar;
1327     real             *shiftvec,*fshift,*x,*f;
1328     real             *fjptrA,*fjptrB,*fjptrC,*fjptrD;
1329     real             scratch[4*DIM];
1330     __m128           fscal,rcutoff,rcutoff2,jidxall;
1331     int              vdwioffset0;
1332     __m128           ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
1333     int              vdwioffset1;
1334     __m128           ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
1335     int              vdwioffset2;
1336     __m128           ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
1337     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
1338     __m128           jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
1339     int              vdwjidx1A,vdwjidx1B,vdwjidx1C,vdwjidx1D;
1340     __m128           jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
1341     int              vdwjidx2A,vdwjidx2B,vdwjidx2C,vdwjidx2D;
1342     __m128           jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
1343     __m128           dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
1344     __m128           dx01,dy01,dz01,rsq01,rinv01,rinvsq01,r01,qq01,c6_01,c12_01;
1345     __m128           dx02,dy02,dz02,rsq02,rinv02,rinvsq02,r02,qq02,c6_02,c12_02;
1346     __m128           dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
1347     __m128           dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11;
1348     __m128           dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12;
1349     __m128           dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
1350     __m128           dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21;
1351     __m128           dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22;
1352     __m128           velec,felec,velecsum,facel,crf,krf,krf2;
1353     real             *charge;
1354     int              nvdwtype;
1355     __m128           rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
1356     int              *vdwtype;
1357     real             *vdwparam;
1358     __m128           one_sixth   = _mm_set1_ps(1.0/6.0);
1359     __m128           one_twelfth = _mm_set1_ps(1.0/12.0);
1360     __m128           c6grid_00;
1361     __m128           c6grid_01;
1362     __m128           c6grid_02;
1363     __m128           c6grid_10;
1364     __m128           c6grid_11;
1365     __m128           c6grid_12;
1366     __m128           c6grid_20;
1367     __m128           c6grid_21;
1368     __m128           c6grid_22;
1369     real             *vdwgridparam;
1370     __m128           ewclj,ewclj2,ewclj6,ewcljrsq,poly,exponent,f6A,f6B,sh_lj_ewald;
1371     __m128           one_half = _mm_set1_ps(0.5);
1372     __m128           minus_one = _mm_set1_ps(-1.0);
1373     __m128i          ewitab;
1374     __m128           ewtabscale,eweps,twoeweps,sh_ewald,ewrt,ewtabhalfspace,ewtabF,ewtabFn,ewtabD,ewtabV;
1375     __m128           beta,beta2,beta3,zeta2,pmecorrF,pmecorrV,rinv3;
1376     real             *ewtab;
1377     __m128           dummy_mask,cutoff_mask;
1378     __m128           signbit = _mm_castsi128_ps( _mm_set1_epi32(0x80000000) );
1379     __m128           one     = _mm_set1_ps(1.0);
1380     __m128           two     = _mm_set1_ps(2.0);
1381     x                = xx[0];
1382     f                = ff[0];
1383
1384     nri              = nlist->nri;
1385     iinr             = nlist->iinr;
1386     jindex           = nlist->jindex;
1387     jjnr             = nlist->jjnr;
1388     shiftidx         = nlist->shift;
1389     gid              = nlist->gid;
1390     shiftvec         = fr->shift_vec[0];
1391     fshift           = fr->fshift[0];
1392     facel            = _mm_set1_ps(fr->ic->epsfac);
1393     charge           = mdatoms->chargeA;
1394     nvdwtype         = fr->ntype;
1395     vdwparam         = fr->nbfp;
1396     vdwtype          = mdatoms->typeA;
1397     vdwgridparam     = fr->ljpme_c6grid;
1398     sh_lj_ewald      = _mm_set1_ps(fr->ic->sh_lj_ewald);
1399     ewclj            = _mm_set1_ps(fr->ic->ewaldcoeff_lj);
1400     ewclj2           = _mm_mul_ps(minus_one,_mm_mul_ps(ewclj,ewclj));
1401
1402     sh_ewald         = _mm_set1_ps(fr->ic->sh_ewald);
1403     beta             = _mm_set1_ps(fr->ic->ewaldcoeff_q);
1404     beta2            = _mm_mul_ps(beta,beta);
1405     beta3            = _mm_mul_ps(beta,beta2);
1406     ewtab            = fr->ic->tabq_coul_F;
1407     ewtabscale       = _mm_set1_ps(fr->ic->tabq_scale);
1408     ewtabhalfspace   = _mm_set1_ps(0.5/fr->ic->tabq_scale);
1409
1410     /* Setup water-specific parameters */
1411     inr              = nlist->iinr[0];
1412     iq0              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+0]));
1413     iq1              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+1]));
1414     iq2              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+2]));
1415     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
1416
1417     jq0              = _mm_set1_ps(charge[inr+0]);
1418     jq1              = _mm_set1_ps(charge[inr+1]);
1419     jq2              = _mm_set1_ps(charge[inr+2]);
1420     vdwjidx0A        = 2*vdwtype[inr+0];
1421     qq00             = _mm_mul_ps(iq0,jq0);
1422     c6_00            = _mm_set1_ps(vdwparam[vdwioffset0+vdwjidx0A]);
1423     c12_00           = _mm_set1_ps(vdwparam[vdwioffset0+vdwjidx0A+1]);
1424     c6grid_00        = _mm_set1_ps(vdwgridparam[vdwioffset0+vdwjidx0A]);
1425     qq01             = _mm_mul_ps(iq0,jq1);
1426     qq02             = _mm_mul_ps(iq0,jq2);
1427     qq10             = _mm_mul_ps(iq1,jq0);
1428     qq11             = _mm_mul_ps(iq1,jq1);
1429     qq12             = _mm_mul_ps(iq1,jq2);
1430     qq20             = _mm_mul_ps(iq2,jq0);
1431     qq21             = _mm_mul_ps(iq2,jq1);
1432     qq22             = _mm_mul_ps(iq2,jq2);
1433
1434     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
1435     rcutoff_scalar   = fr->ic->rcoulomb;
1436     rcutoff          = _mm_set1_ps(rcutoff_scalar);
1437     rcutoff2         = _mm_mul_ps(rcutoff,rcutoff);
1438
1439     sh_vdw_invrcut6  = _mm_set1_ps(fr->ic->sh_invrc6);
1440     rvdw             = _mm_set1_ps(fr->ic->rvdw);
1441
1442     /* Avoid stupid compiler warnings */
1443     jnrA = jnrB = jnrC = jnrD = 0;
1444     j_coord_offsetA = 0;
1445     j_coord_offsetB = 0;
1446     j_coord_offsetC = 0;
1447     j_coord_offsetD = 0;
1448
1449     outeriter        = 0;
1450     inneriter        = 0;
1451
1452     for(iidx=0;iidx<4*DIM;iidx++)
1453     {
1454         scratch[iidx] = 0.0;
1455     }
1456
1457     /* Start outer loop over neighborlists */
1458     for(iidx=0; iidx<nri; iidx++)
1459     {
1460         /* Load shift vector for this list */
1461         i_shift_offset   = DIM*shiftidx[iidx];
1462
1463         /* Load limits for loop over neighbors */
1464         j_index_start    = jindex[iidx];
1465         j_index_end      = jindex[iidx+1];
1466
1467         /* Get outer coordinate index */
1468         inr              = iinr[iidx];
1469         i_coord_offset   = DIM*inr;
1470
1471         /* Load i particle coords and add shift vector */
1472         gmx_mm_load_shift_and_3rvec_broadcast_ps(shiftvec+i_shift_offset,x+i_coord_offset,
1473                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2);
1474
1475         fix0             = _mm_setzero_ps();
1476         fiy0             = _mm_setzero_ps();
1477         fiz0             = _mm_setzero_ps();
1478         fix1             = _mm_setzero_ps();
1479         fiy1             = _mm_setzero_ps();
1480         fiz1             = _mm_setzero_ps();
1481         fix2             = _mm_setzero_ps();
1482         fiy2             = _mm_setzero_ps();
1483         fiz2             = _mm_setzero_ps();
1484
1485         /* Start inner kernel loop */
1486         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
1487         {
1488
1489             /* Get j neighbor index, and coordinate index */
1490             jnrA             = jjnr[jidx];
1491             jnrB             = jjnr[jidx+1];
1492             jnrC             = jjnr[jidx+2];
1493             jnrD             = jjnr[jidx+3];
1494             j_coord_offsetA  = DIM*jnrA;
1495             j_coord_offsetB  = DIM*jnrB;
1496             j_coord_offsetC  = DIM*jnrC;
1497             j_coord_offsetD  = DIM*jnrD;
1498
1499             /* load j atom coordinates */
1500             gmx_mm_load_3rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
1501                                               x+j_coord_offsetC,x+j_coord_offsetD,
1502                                               &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,&jy2,&jz2);
1503
1504             /* Calculate displacement vector */
1505             dx00             = _mm_sub_ps(ix0,jx0);
1506             dy00             = _mm_sub_ps(iy0,jy0);
1507             dz00             = _mm_sub_ps(iz0,jz0);
1508             dx01             = _mm_sub_ps(ix0,jx1);
1509             dy01             = _mm_sub_ps(iy0,jy1);
1510             dz01             = _mm_sub_ps(iz0,jz1);
1511             dx02             = _mm_sub_ps(ix0,jx2);
1512             dy02             = _mm_sub_ps(iy0,jy2);
1513             dz02             = _mm_sub_ps(iz0,jz2);
1514             dx10             = _mm_sub_ps(ix1,jx0);
1515             dy10             = _mm_sub_ps(iy1,jy0);
1516             dz10             = _mm_sub_ps(iz1,jz0);
1517             dx11             = _mm_sub_ps(ix1,jx1);
1518             dy11             = _mm_sub_ps(iy1,jy1);
1519             dz11             = _mm_sub_ps(iz1,jz1);
1520             dx12             = _mm_sub_ps(ix1,jx2);
1521             dy12             = _mm_sub_ps(iy1,jy2);
1522             dz12             = _mm_sub_ps(iz1,jz2);
1523             dx20             = _mm_sub_ps(ix2,jx0);
1524             dy20             = _mm_sub_ps(iy2,jy0);
1525             dz20             = _mm_sub_ps(iz2,jz0);
1526             dx21             = _mm_sub_ps(ix2,jx1);
1527             dy21             = _mm_sub_ps(iy2,jy1);
1528             dz21             = _mm_sub_ps(iz2,jz1);
1529             dx22             = _mm_sub_ps(ix2,jx2);
1530             dy22             = _mm_sub_ps(iy2,jy2);
1531             dz22             = _mm_sub_ps(iz2,jz2);
1532
1533             /* Calculate squared distance and things based on it */
1534             rsq00            = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
1535             rsq01            = gmx_mm_calc_rsq_ps(dx01,dy01,dz01);
1536             rsq02            = gmx_mm_calc_rsq_ps(dx02,dy02,dz02);
1537             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
1538             rsq11            = gmx_mm_calc_rsq_ps(dx11,dy11,dz11);
1539             rsq12            = gmx_mm_calc_rsq_ps(dx12,dy12,dz12);
1540             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
1541             rsq21            = gmx_mm_calc_rsq_ps(dx21,dy21,dz21);
1542             rsq22            = gmx_mm_calc_rsq_ps(dx22,dy22,dz22);
1543
1544             rinv00           = avx128fma_invsqrt_f(rsq00);
1545             rinv01           = avx128fma_invsqrt_f(rsq01);
1546             rinv02           = avx128fma_invsqrt_f(rsq02);
1547             rinv10           = avx128fma_invsqrt_f(rsq10);
1548             rinv11           = avx128fma_invsqrt_f(rsq11);
1549             rinv12           = avx128fma_invsqrt_f(rsq12);
1550             rinv20           = avx128fma_invsqrt_f(rsq20);
1551             rinv21           = avx128fma_invsqrt_f(rsq21);
1552             rinv22           = avx128fma_invsqrt_f(rsq22);
1553
1554             rinvsq00         = _mm_mul_ps(rinv00,rinv00);
1555             rinvsq01         = _mm_mul_ps(rinv01,rinv01);
1556             rinvsq02         = _mm_mul_ps(rinv02,rinv02);
1557             rinvsq10         = _mm_mul_ps(rinv10,rinv10);
1558             rinvsq11         = _mm_mul_ps(rinv11,rinv11);
1559             rinvsq12         = _mm_mul_ps(rinv12,rinv12);
1560             rinvsq20         = _mm_mul_ps(rinv20,rinv20);
1561             rinvsq21         = _mm_mul_ps(rinv21,rinv21);
1562             rinvsq22         = _mm_mul_ps(rinv22,rinv22);
1563
1564             fjx0             = _mm_setzero_ps();
1565             fjy0             = _mm_setzero_ps();
1566             fjz0             = _mm_setzero_ps();
1567             fjx1             = _mm_setzero_ps();
1568             fjy1             = _mm_setzero_ps();
1569             fjz1             = _mm_setzero_ps();
1570             fjx2             = _mm_setzero_ps();
1571             fjy2             = _mm_setzero_ps();
1572             fjz2             = _mm_setzero_ps();
1573
1574             /**************************
1575              * CALCULATE INTERACTIONS *
1576              **************************/
1577
1578             if (gmx_mm_any_lt(rsq00,rcutoff2))
1579             {
1580
1581             r00              = _mm_mul_ps(rsq00,rinv00);
1582
1583             /* EWALD ELECTROSTATICS */
1584
1585             /* Analytical PME correction */
1586             zeta2            = _mm_mul_ps(beta2,rsq00);
1587             rinv3            = _mm_mul_ps(rinvsq00,rinv00);
1588             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
1589             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
1590             felec            = _mm_mul_ps(qq00,felec);
1591
1592             /* Analytical LJ-PME */
1593             rinvsix          = _mm_mul_ps(_mm_mul_ps(rinvsq00,rinvsq00),rinvsq00);
1594             ewcljrsq         = _mm_mul_ps(ewclj2,rsq00);
1595             ewclj6           = _mm_mul_ps(ewclj2,_mm_mul_ps(ewclj2,ewclj2));
1596             exponent         = avx128fma_exp_f(ewcljrsq);
1597             /* poly = exp(-(beta*r)^2) * (1 + (beta*r)^2 + (beta*r)^4 /2) */
1598             poly             = _mm_mul_ps(exponent,_mm_macc_ps(_mm_mul_ps(ewcljrsq,ewcljrsq),one_half,_mm_sub_ps(one,ewcljrsq)));
1599             /* f6A = 6 * C6grid * (1 - poly) */
1600             f6A              = _mm_mul_ps(c6grid_00,_mm_sub_ps(one,poly));
1601             /* f6B = C6grid * exponent * beta^6 */
1602             f6B              = _mm_mul_ps(_mm_mul_ps(c6grid_00,one_sixth),_mm_mul_ps(exponent,ewclj6));
1603             /* fvdw = 12*C12/r13 - ((6*C6 - f6A)/r6 + f6B)/r */
1604             fvdw              = _mm_mul_ps(_mm_macc_ps(_mm_msub_ps(c12_00,rinvsix,_mm_sub_ps(c6_00,f6A)),rinvsix,f6B),rinvsq00);
1605
1606             cutoff_mask      = _mm_cmplt_ps(rsq00,rcutoff2);
1607
1608             fscal            = _mm_add_ps(felec,fvdw);
1609
1610             fscal            = _mm_and_ps(fscal,cutoff_mask);
1611
1612              /* Update vectorial force */
1613             fix0             = _mm_macc_ps(dx00,fscal,fix0);
1614             fiy0             = _mm_macc_ps(dy00,fscal,fiy0);
1615             fiz0             = _mm_macc_ps(dz00,fscal,fiz0);
1616
1617             fjx0             = _mm_macc_ps(dx00,fscal,fjx0);
1618             fjy0             = _mm_macc_ps(dy00,fscal,fjy0);
1619             fjz0             = _mm_macc_ps(dz00,fscal,fjz0);
1620
1621             }
1622
1623             /**************************
1624              * CALCULATE INTERACTIONS *
1625              **************************/
1626
1627             if (gmx_mm_any_lt(rsq01,rcutoff2))
1628             {
1629
1630             r01              = _mm_mul_ps(rsq01,rinv01);
1631
1632             /* EWALD ELECTROSTATICS */
1633
1634             /* Analytical PME correction */
1635             zeta2            = _mm_mul_ps(beta2,rsq01);
1636             rinv3            = _mm_mul_ps(rinvsq01,rinv01);
1637             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
1638             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
1639             felec            = _mm_mul_ps(qq01,felec);
1640
1641             cutoff_mask      = _mm_cmplt_ps(rsq01,rcutoff2);
1642
1643             fscal            = felec;
1644
1645             fscal            = _mm_and_ps(fscal,cutoff_mask);
1646
1647              /* Update vectorial force */
1648             fix0             = _mm_macc_ps(dx01,fscal,fix0);
1649             fiy0             = _mm_macc_ps(dy01,fscal,fiy0);
1650             fiz0             = _mm_macc_ps(dz01,fscal,fiz0);
1651
1652             fjx1             = _mm_macc_ps(dx01,fscal,fjx1);
1653             fjy1             = _mm_macc_ps(dy01,fscal,fjy1);
1654             fjz1             = _mm_macc_ps(dz01,fscal,fjz1);
1655
1656             }
1657
1658             /**************************
1659              * CALCULATE INTERACTIONS *
1660              **************************/
1661
1662             if (gmx_mm_any_lt(rsq02,rcutoff2))
1663             {
1664
1665             r02              = _mm_mul_ps(rsq02,rinv02);
1666
1667             /* EWALD ELECTROSTATICS */
1668
1669             /* Analytical PME correction */
1670             zeta2            = _mm_mul_ps(beta2,rsq02);
1671             rinv3            = _mm_mul_ps(rinvsq02,rinv02);
1672             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
1673             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
1674             felec            = _mm_mul_ps(qq02,felec);
1675
1676             cutoff_mask      = _mm_cmplt_ps(rsq02,rcutoff2);
1677
1678             fscal            = felec;
1679
1680             fscal            = _mm_and_ps(fscal,cutoff_mask);
1681
1682              /* Update vectorial force */
1683             fix0             = _mm_macc_ps(dx02,fscal,fix0);
1684             fiy0             = _mm_macc_ps(dy02,fscal,fiy0);
1685             fiz0             = _mm_macc_ps(dz02,fscal,fiz0);
1686
1687             fjx2             = _mm_macc_ps(dx02,fscal,fjx2);
1688             fjy2             = _mm_macc_ps(dy02,fscal,fjy2);
1689             fjz2             = _mm_macc_ps(dz02,fscal,fjz2);
1690
1691             }
1692
1693             /**************************
1694              * CALCULATE INTERACTIONS *
1695              **************************/
1696
1697             if (gmx_mm_any_lt(rsq10,rcutoff2))
1698             {
1699
1700             r10              = _mm_mul_ps(rsq10,rinv10);
1701
1702             /* EWALD ELECTROSTATICS */
1703
1704             /* Analytical PME correction */
1705             zeta2            = _mm_mul_ps(beta2,rsq10);
1706             rinv3            = _mm_mul_ps(rinvsq10,rinv10);
1707             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
1708             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
1709             felec            = _mm_mul_ps(qq10,felec);
1710
1711             cutoff_mask      = _mm_cmplt_ps(rsq10,rcutoff2);
1712
1713             fscal            = felec;
1714
1715             fscal            = _mm_and_ps(fscal,cutoff_mask);
1716
1717              /* Update vectorial force */
1718             fix1             = _mm_macc_ps(dx10,fscal,fix1);
1719             fiy1             = _mm_macc_ps(dy10,fscal,fiy1);
1720             fiz1             = _mm_macc_ps(dz10,fscal,fiz1);
1721
1722             fjx0             = _mm_macc_ps(dx10,fscal,fjx0);
1723             fjy0             = _mm_macc_ps(dy10,fscal,fjy0);
1724             fjz0             = _mm_macc_ps(dz10,fscal,fjz0);
1725
1726             }
1727
1728             /**************************
1729              * CALCULATE INTERACTIONS *
1730              **************************/
1731
1732             if (gmx_mm_any_lt(rsq11,rcutoff2))
1733             {
1734
1735             r11              = _mm_mul_ps(rsq11,rinv11);
1736
1737             /* EWALD ELECTROSTATICS */
1738
1739             /* Analytical PME correction */
1740             zeta2            = _mm_mul_ps(beta2,rsq11);
1741             rinv3            = _mm_mul_ps(rinvsq11,rinv11);
1742             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
1743             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
1744             felec            = _mm_mul_ps(qq11,felec);
1745
1746             cutoff_mask      = _mm_cmplt_ps(rsq11,rcutoff2);
1747
1748             fscal            = felec;
1749
1750             fscal            = _mm_and_ps(fscal,cutoff_mask);
1751
1752              /* Update vectorial force */
1753             fix1             = _mm_macc_ps(dx11,fscal,fix1);
1754             fiy1             = _mm_macc_ps(dy11,fscal,fiy1);
1755             fiz1             = _mm_macc_ps(dz11,fscal,fiz1);
1756
1757             fjx1             = _mm_macc_ps(dx11,fscal,fjx1);
1758             fjy1             = _mm_macc_ps(dy11,fscal,fjy1);
1759             fjz1             = _mm_macc_ps(dz11,fscal,fjz1);
1760
1761             }
1762
1763             /**************************
1764              * CALCULATE INTERACTIONS *
1765              **************************/
1766
1767             if (gmx_mm_any_lt(rsq12,rcutoff2))
1768             {
1769
1770             r12              = _mm_mul_ps(rsq12,rinv12);
1771
1772             /* EWALD ELECTROSTATICS */
1773
1774             /* Analytical PME correction */
1775             zeta2            = _mm_mul_ps(beta2,rsq12);
1776             rinv3            = _mm_mul_ps(rinvsq12,rinv12);
1777             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
1778             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
1779             felec            = _mm_mul_ps(qq12,felec);
1780
1781             cutoff_mask      = _mm_cmplt_ps(rsq12,rcutoff2);
1782
1783             fscal            = felec;
1784
1785             fscal            = _mm_and_ps(fscal,cutoff_mask);
1786
1787              /* Update vectorial force */
1788             fix1             = _mm_macc_ps(dx12,fscal,fix1);
1789             fiy1             = _mm_macc_ps(dy12,fscal,fiy1);
1790             fiz1             = _mm_macc_ps(dz12,fscal,fiz1);
1791
1792             fjx2             = _mm_macc_ps(dx12,fscal,fjx2);
1793             fjy2             = _mm_macc_ps(dy12,fscal,fjy2);
1794             fjz2             = _mm_macc_ps(dz12,fscal,fjz2);
1795
1796             }
1797
1798             /**************************
1799              * CALCULATE INTERACTIONS *
1800              **************************/
1801
1802             if (gmx_mm_any_lt(rsq20,rcutoff2))
1803             {
1804
1805             r20              = _mm_mul_ps(rsq20,rinv20);
1806
1807             /* EWALD ELECTROSTATICS */
1808
1809             /* Analytical PME correction */
1810             zeta2            = _mm_mul_ps(beta2,rsq20);
1811             rinv3            = _mm_mul_ps(rinvsq20,rinv20);
1812             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
1813             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
1814             felec            = _mm_mul_ps(qq20,felec);
1815
1816             cutoff_mask      = _mm_cmplt_ps(rsq20,rcutoff2);
1817
1818             fscal            = felec;
1819
1820             fscal            = _mm_and_ps(fscal,cutoff_mask);
1821
1822              /* Update vectorial force */
1823             fix2             = _mm_macc_ps(dx20,fscal,fix2);
1824             fiy2             = _mm_macc_ps(dy20,fscal,fiy2);
1825             fiz2             = _mm_macc_ps(dz20,fscal,fiz2);
1826
1827             fjx0             = _mm_macc_ps(dx20,fscal,fjx0);
1828             fjy0             = _mm_macc_ps(dy20,fscal,fjy0);
1829             fjz0             = _mm_macc_ps(dz20,fscal,fjz0);
1830
1831             }
1832
1833             /**************************
1834              * CALCULATE INTERACTIONS *
1835              **************************/
1836
1837             if (gmx_mm_any_lt(rsq21,rcutoff2))
1838             {
1839
1840             r21              = _mm_mul_ps(rsq21,rinv21);
1841
1842             /* EWALD ELECTROSTATICS */
1843
1844             /* Analytical PME correction */
1845             zeta2            = _mm_mul_ps(beta2,rsq21);
1846             rinv3            = _mm_mul_ps(rinvsq21,rinv21);
1847             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
1848             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
1849             felec            = _mm_mul_ps(qq21,felec);
1850
1851             cutoff_mask      = _mm_cmplt_ps(rsq21,rcutoff2);
1852
1853             fscal            = felec;
1854
1855             fscal            = _mm_and_ps(fscal,cutoff_mask);
1856
1857              /* Update vectorial force */
1858             fix2             = _mm_macc_ps(dx21,fscal,fix2);
1859             fiy2             = _mm_macc_ps(dy21,fscal,fiy2);
1860             fiz2             = _mm_macc_ps(dz21,fscal,fiz2);
1861
1862             fjx1             = _mm_macc_ps(dx21,fscal,fjx1);
1863             fjy1             = _mm_macc_ps(dy21,fscal,fjy1);
1864             fjz1             = _mm_macc_ps(dz21,fscal,fjz1);
1865
1866             }
1867
1868             /**************************
1869              * CALCULATE INTERACTIONS *
1870              **************************/
1871
1872             if (gmx_mm_any_lt(rsq22,rcutoff2))
1873             {
1874
1875             r22              = _mm_mul_ps(rsq22,rinv22);
1876
1877             /* EWALD ELECTROSTATICS */
1878
1879             /* Analytical PME correction */
1880             zeta2            = _mm_mul_ps(beta2,rsq22);
1881             rinv3            = _mm_mul_ps(rinvsq22,rinv22);
1882             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
1883             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
1884             felec            = _mm_mul_ps(qq22,felec);
1885
1886             cutoff_mask      = _mm_cmplt_ps(rsq22,rcutoff2);
1887
1888             fscal            = felec;
1889
1890             fscal            = _mm_and_ps(fscal,cutoff_mask);
1891
1892              /* Update vectorial force */
1893             fix2             = _mm_macc_ps(dx22,fscal,fix2);
1894             fiy2             = _mm_macc_ps(dy22,fscal,fiy2);
1895             fiz2             = _mm_macc_ps(dz22,fscal,fiz2);
1896
1897             fjx2             = _mm_macc_ps(dx22,fscal,fjx2);
1898             fjy2             = _mm_macc_ps(dy22,fscal,fjy2);
1899             fjz2             = _mm_macc_ps(dz22,fscal,fjz2);
1900
1901             }
1902
1903             fjptrA             = f+j_coord_offsetA;
1904             fjptrB             = f+j_coord_offsetB;
1905             fjptrC             = f+j_coord_offsetC;
1906             fjptrD             = f+j_coord_offsetD;
1907
1908             gmx_mm_decrement_3rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,
1909                                                    fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,fjx2,fjy2,fjz2);
1910
1911             /* Inner loop uses 300 flops */
1912         }
1913
1914         if(jidx<j_index_end)
1915         {
1916
1917             /* Get j neighbor index, and coordinate index */
1918             jnrlistA         = jjnr[jidx];
1919             jnrlistB         = jjnr[jidx+1];
1920             jnrlistC         = jjnr[jidx+2];
1921             jnrlistD         = jjnr[jidx+3];
1922             /* Sign of each element will be negative for non-real atoms.
1923              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
1924              * so use it as val = _mm_andnot_ps(mask,val) to clear dummy entries.
1925              */
1926             dummy_mask = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
1927             jnrA       = (jnrlistA>=0) ? jnrlistA : 0;
1928             jnrB       = (jnrlistB>=0) ? jnrlistB : 0;
1929             jnrC       = (jnrlistC>=0) ? jnrlistC : 0;
1930             jnrD       = (jnrlistD>=0) ? jnrlistD : 0;
1931             j_coord_offsetA  = DIM*jnrA;
1932             j_coord_offsetB  = DIM*jnrB;
1933             j_coord_offsetC  = DIM*jnrC;
1934             j_coord_offsetD  = DIM*jnrD;
1935
1936             /* load j atom coordinates */
1937             gmx_mm_load_3rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
1938                                               x+j_coord_offsetC,x+j_coord_offsetD,
1939                                               &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,&jy2,&jz2);
1940
1941             /* Calculate displacement vector */
1942             dx00             = _mm_sub_ps(ix0,jx0);
1943             dy00             = _mm_sub_ps(iy0,jy0);
1944             dz00             = _mm_sub_ps(iz0,jz0);
1945             dx01             = _mm_sub_ps(ix0,jx1);
1946             dy01             = _mm_sub_ps(iy0,jy1);
1947             dz01             = _mm_sub_ps(iz0,jz1);
1948             dx02             = _mm_sub_ps(ix0,jx2);
1949             dy02             = _mm_sub_ps(iy0,jy2);
1950             dz02             = _mm_sub_ps(iz0,jz2);
1951             dx10             = _mm_sub_ps(ix1,jx0);
1952             dy10             = _mm_sub_ps(iy1,jy0);
1953             dz10             = _mm_sub_ps(iz1,jz0);
1954             dx11             = _mm_sub_ps(ix1,jx1);
1955             dy11             = _mm_sub_ps(iy1,jy1);
1956             dz11             = _mm_sub_ps(iz1,jz1);
1957             dx12             = _mm_sub_ps(ix1,jx2);
1958             dy12             = _mm_sub_ps(iy1,jy2);
1959             dz12             = _mm_sub_ps(iz1,jz2);
1960             dx20             = _mm_sub_ps(ix2,jx0);
1961             dy20             = _mm_sub_ps(iy2,jy0);
1962             dz20             = _mm_sub_ps(iz2,jz0);
1963             dx21             = _mm_sub_ps(ix2,jx1);
1964             dy21             = _mm_sub_ps(iy2,jy1);
1965             dz21             = _mm_sub_ps(iz2,jz1);
1966             dx22             = _mm_sub_ps(ix2,jx2);
1967             dy22             = _mm_sub_ps(iy2,jy2);
1968             dz22             = _mm_sub_ps(iz2,jz2);
1969
1970             /* Calculate squared distance and things based on it */
1971             rsq00            = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
1972             rsq01            = gmx_mm_calc_rsq_ps(dx01,dy01,dz01);
1973             rsq02            = gmx_mm_calc_rsq_ps(dx02,dy02,dz02);
1974             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
1975             rsq11            = gmx_mm_calc_rsq_ps(dx11,dy11,dz11);
1976             rsq12            = gmx_mm_calc_rsq_ps(dx12,dy12,dz12);
1977             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
1978             rsq21            = gmx_mm_calc_rsq_ps(dx21,dy21,dz21);
1979             rsq22            = gmx_mm_calc_rsq_ps(dx22,dy22,dz22);
1980
1981             rinv00           = avx128fma_invsqrt_f(rsq00);
1982             rinv01           = avx128fma_invsqrt_f(rsq01);
1983             rinv02           = avx128fma_invsqrt_f(rsq02);
1984             rinv10           = avx128fma_invsqrt_f(rsq10);
1985             rinv11           = avx128fma_invsqrt_f(rsq11);
1986             rinv12           = avx128fma_invsqrt_f(rsq12);
1987             rinv20           = avx128fma_invsqrt_f(rsq20);
1988             rinv21           = avx128fma_invsqrt_f(rsq21);
1989             rinv22           = avx128fma_invsqrt_f(rsq22);
1990
1991             rinvsq00         = _mm_mul_ps(rinv00,rinv00);
1992             rinvsq01         = _mm_mul_ps(rinv01,rinv01);
1993             rinvsq02         = _mm_mul_ps(rinv02,rinv02);
1994             rinvsq10         = _mm_mul_ps(rinv10,rinv10);
1995             rinvsq11         = _mm_mul_ps(rinv11,rinv11);
1996             rinvsq12         = _mm_mul_ps(rinv12,rinv12);
1997             rinvsq20         = _mm_mul_ps(rinv20,rinv20);
1998             rinvsq21         = _mm_mul_ps(rinv21,rinv21);
1999             rinvsq22         = _mm_mul_ps(rinv22,rinv22);
2000
2001             fjx0             = _mm_setzero_ps();
2002             fjy0             = _mm_setzero_ps();
2003             fjz0             = _mm_setzero_ps();
2004             fjx1             = _mm_setzero_ps();
2005             fjy1             = _mm_setzero_ps();
2006             fjz1             = _mm_setzero_ps();
2007             fjx2             = _mm_setzero_ps();
2008             fjy2             = _mm_setzero_ps();
2009             fjz2             = _mm_setzero_ps();
2010
2011             /**************************
2012              * CALCULATE INTERACTIONS *
2013              **************************/
2014
2015             if (gmx_mm_any_lt(rsq00,rcutoff2))
2016             {
2017
2018             r00              = _mm_mul_ps(rsq00,rinv00);
2019             r00              = _mm_andnot_ps(dummy_mask,r00);
2020
2021             /* EWALD ELECTROSTATICS */
2022
2023             /* Analytical PME correction */
2024             zeta2            = _mm_mul_ps(beta2,rsq00);
2025             rinv3            = _mm_mul_ps(rinvsq00,rinv00);
2026             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
2027             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
2028             felec            = _mm_mul_ps(qq00,felec);
2029
2030             /* Analytical LJ-PME */
2031             rinvsix          = _mm_mul_ps(_mm_mul_ps(rinvsq00,rinvsq00),rinvsq00);
2032             ewcljrsq         = _mm_mul_ps(ewclj2,rsq00);
2033             ewclj6           = _mm_mul_ps(ewclj2,_mm_mul_ps(ewclj2,ewclj2));
2034             exponent         = avx128fma_exp_f(ewcljrsq);
2035             /* poly = exp(-(beta*r)^2) * (1 + (beta*r)^2 + (beta*r)^4 /2) */
2036             poly             = _mm_mul_ps(exponent,_mm_macc_ps(_mm_mul_ps(ewcljrsq,ewcljrsq),one_half,_mm_sub_ps(one,ewcljrsq)));
2037             /* f6A = 6 * C6grid * (1 - poly) */
2038             f6A              = _mm_mul_ps(c6grid_00,_mm_sub_ps(one,poly));
2039             /* f6B = C6grid * exponent * beta^6 */
2040             f6B              = _mm_mul_ps(_mm_mul_ps(c6grid_00,one_sixth),_mm_mul_ps(exponent,ewclj6));
2041             /* fvdw = 12*C12/r13 - ((6*C6 - f6A)/r6 + f6B)/r */
2042             fvdw              = _mm_mul_ps(_mm_macc_ps(_mm_msub_ps(c12_00,rinvsix,_mm_sub_ps(c6_00,f6A)),rinvsix,f6B),rinvsq00);
2043
2044             cutoff_mask      = _mm_cmplt_ps(rsq00,rcutoff2);
2045
2046             fscal            = _mm_add_ps(felec,fvdw);
2047
2048             fscal            = _mm_and_ps(fscal,cutoff_mask);
2049
2050             fscal            = _mm_andnot_ps(dummy_mask,fscal);
2051
2052              /* Update vectorial force */
2053             fix0             = _mm_macc_ps(dx00,fscal,fix0);
2054             fiy0             = _mm_macc_ps(dy00,fscal,fiy0);
2055             fiz0             = _mm_macc_ps(dz00,fscal,fiz0);
2056
2057             fjx0             = _mm_macc_ps(dx00,fscal,fjx0);
2058             fjy0             = _mm_macc_ps(dy00,fscal,fjy0);
2059             fjz0             = _mm_macc_ps(dz00,fscal,fjz0);
2060
2061             }
2062
2063             /**************************
2064              * CALCULATE INTERACTIONS *
2065              **************************/
2066
2067             if (gmx_mm_any_lt(rsq01,rcutoff2))
2068             {
2069
2070             r01              = _mm_mul_ps(rsq01,rinv01);
2071             r01              = _mm_andnot_ps(dummy_mask,r01);
2072
2073             /* EWALD ELECTROSTATICS */
2074
2075             /* Analytical PME correction */
2076             zeta2            = _mm_mul_ps(beta2,rsq01);
2077             rinv3            = _mm_mul_ps(rinvsq01,rinv01);
2078             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
2079             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
2080             felec            = _mm_mul_ps(qq01,felec);
2081
2082             cutoff_mask      = _mm_cmplt_ps(rsq01,rcutoff2);
2083
2084             fscal            = felec;
2085
2086             fscal            = _mm_and_ps(fscal,cutoff_mask);
2087
2088             fscal            = _mm_andnot_ps(dummy_mask,fscal);
2089
2090              /* Update vectorial force */
2091             fix0             = _mm_macc_ps(dx01,fscal,fix0);
2092             fiy0             = _mm_macc_ps(dy01,fscal,fiy0);
2093             fiz0             = _mm_macc_ps(dz01,fscal,fiz0);
2094
2095             fjx1             = _mm_macc_ps(dx01,fscal,fjx1);
2096             fjy1             = _mm_macc_ps(dy01,fscal,fjy1);
2097             fjz1             = _mm_macc_ps(dz01,fscal,fjz1);
2098
2099             }
2100
2101             /**************************
2102              * CALCULATE INTERACTIONS *
2103              **************************/
2104
2105             if (gmx_mm_any_lt(rsq02,rcutoff2))
2106             {
2107
2108             r02              = _mm_mul_ps(rsq02,rinv02);
2109             r02              = _mm_andnot_ps(dummy_mask,r02);
2110
2111             /* EWALD ELECTROSTATICS */
2112
2113             /* Analytical PME correction */
2114             zeta2            = _mm_mul_ps(beta2,rsq02);
2115             rinv3            = _mm_mul_ps(rinvsq02,rinv02);
2116             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
2117             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
2118             felec            = _mm_mul_ps(qq02,felec);
2119
2120             cutoff_mask      = _mm_cmplt_ps(rsq02,rcutoff2);
2121
2122             fscal            = felec;
2123
2124             fscal            = _mm_and_ps(fscal,cutoff_mask);
2125
2126             fscal            = _mm_andnot_ps(dummy_mask,fscal);
2127
2128              /* Update vectorial force */
2129             fix0             = _mm_macc_ps(dx02,fscal,fix0);
2130             fiy0             = _mm_macc_ps(dy02,fscal,fiy0);
2131             fiz0             = _mm_macc_ps(dz02,fscal,fiz0);
2132
2133             fjx2             = _mm_macc_ps(dx02,fscal,fjx2);
2134             fjy2             = _mm_macc_ps(dy02,fscal,fjy2);
2135             fjz2             = _mm_macc_ps(dz02,fscal,fjz2);
2136
2137             }
2138
2139             /**************************
2140              * CALCULATE INTERACTIONS *
2141              **************************/
2142
2143             if (gmx_mm_any_lt(rsq10,rcutoff2))
2144             {
2145
2146             r10              = _mm_mul_ps(rsq10,rinv10);
2147             r10              = _mm_andnot_ps(dummy_mask,r10);
2148
2149             /* EWALD ELECTROSTATICS */
2150
2151             /* Analytical PME correction */
2152             zeta2            = _mm_mul_ps(beta2,rsq10);
2153             rinv3            = _mm_mul_ps(rinvsq10,rinv10);
2154             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
2155             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
2156             felec            = _mm_mul_ps(qq10,felec);
2157
2158             cutoff_mask      = _mm_cmplt_ps(rsq10,rcutoff2);
2159
2160             fscal            = felec;
2161
2162             fscal            = _mm_and_ps(fscal,cutoff_mask);
2163
2164             fscal            = _mm_andnot_ps(dummy_mask,fscal);
2165
2166              /* Update vectorial force */
2167             fix1             = _mm_macc_ps(dx10,fscal,fix1);
2168             fiy1             = _mm_macc_ps(dy10,fscal,fiy1);
2169             fiz1             = _mm_macc_ps(dz10,fscal,fiz1);
2170
2171             fjx0             = _mm_macc_ps(dx10,fscal,fjx0);
2172             fjy0             = _mm_macc_ps(dy10,fscal,fjy0);
2173             fjz0             = _mm_macc_ps(dz10,fscal,fjz0);
2174
2175             }
2176
2177             /**************************
2178              * CALCULATE INTERACTIONS *
2179              **************************/
2180
2181             if (gmx_mm_any_lt(rsq11,rcutoff2))
2182             {
2183
2184             r11              = _mm_mul_ps(rsq11,rinv11);
2185             r11              = _mm_andnot_ps(dummy_mask,r11);
2186
2187             /* EWALD ELECTROSTATICS */
2188
2189             /* Analytical PME correction */
2190             zeta2            = _mm_mul_ps(beta2,rsq11);
2191             rinv3            = _mm_mul_ps(rinvsq11,rinv11);
2192             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
2193             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
2194             felec            = _mm_mul_ps(qq11,felec);
2195
2196             cutoff_mask      = _mm_cmplt_ps(rsq11,rcutoff2);
2197
2198             fscal            = felec;
2199
2200             fscal            = _mm_and_ps(fscal,cutoff_mask);
2201
2202             fscal            = _mm_andnot_ps(dummy_mask,fscal);
2203
2204              /* Update vectorial force */
2205             fix1             = _mm_macc_ps(dx11,fscal,fix1);
2206             fiy1             = _mm_macc_ps(dy11,fscal,fiy1);
2207             fiz1             = _mm_macc_ps(dz11,fscal,fiz1);
2208
2209             fjx1             = _mm_macc_ps(dx11,fscal,fjx1);
2210             fjy1             = _mm_macc_ps(dy11,fscal,fjy1);
2211             fjz1             = _mm_macc_ps(dz11,fscal,fjz1);
2212
2213             }
2214
2215             /**************************
2216              * CALCULATE INTERACTIONS *
2217              **************************/
2218
2219             if (gmx_mm_any_lt(rsq12,rcutoff2))
2220             {
2221
2222             r12              = _mm_mul_ps(rsq12,rinv12);
2223             r12              = _mm_andnot_ps(dummy_mask,r12);
2224
2225             /* EWALD ELECTROSTATICS */
2226
2227             /* Analytical PME correction */
2228             zeta2            = _mm_mul_ps(beta2,rsq12);
2229             rinv3            = _mm_mul_ps(rinvsq12,rinv12);
2230             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
2231             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
2232             felec            = _mm_mul_ps(qq12,felec);
2233
2234             cutoff_mask      = _mm_cmplt_ps(rsq12,rcutoff2);
2235
2236             fscal            = felec;
2237
2238             fscal            = _mm_and_ps(fscal,cutoff_mask);
2239
2240             fscal            = _mm_andnot_ps(dummy_mask,fscal);
2241
2242              /* Update vectorial force */
2243             fix1             = _mm_macc_ps(dx12,fscal,fix1);
2244             fiy1             = _mm_macc_ps(dy12,fscal,fiy1);
2245             fiz1             = _mm_macc_ps(dz12,fscal,fiz1);
2246
2247             fjx2             = _mm_macc_ps(dx12,fscal,fjx2);
2248             fjy2             = _mm_macc_ps(dy12,fscal,fjy2);
2249             fjz2             = _mm_macc_ps(dz12,fscal,fjz2);
2250
2251             }
2252
2253             /**************************
2254              * CALCULATE INTERACTIONS *
2255              **************************/
2256
2257             if (gmx_mm_any_lt(rsq20,rcutoff2))
2258             {
2259
2260             r20              = _mm_mul_ps(rsq20,rinv20);
2261             r20              = _mm_andnot_ps(dummy_mask,r20);
2262
2263             /* EWALD ELECTROSTATICS */
2264
2265             /* Analytical PME correction */
2266             zeta2            = _mm_mul_ps(beta2,rsq20);
2267             rinv3            = _mm_mul_ps(rinvsq20,rinv20);
2268             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
2269             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
2270             felec            = _mm_mul_ps(qq20,felec);
2271
2272             cutoff_mask      = _mm_cmplt_ps(rsq20,rcutoff2);
2273
2274             fscal            = felec;
2275
2276             fscal            = _mm_and_ps(fscal,cutoff_mask);
2277
2278             fscal            = _mm_andnot_ps(dummy_mask,fscal);
2279
2280              /* Update vectorial force */
2281             fix2             = _mm_macc_ps(dx20,fscal,fix2);
2282             fiy2             = _mm_macc_ps(dy20,fscal,fiy2);
2283             fiz2             = _mm_macc_ps(dz20,fscal,fiz2);
2284
2285             fjx0             = _mm_macc_ps(dx20,fscal,fjx0);
2286             fjy0             = _mm_macc_ps(dy20,fscal,fjy0);
2287             fjz0             = _mm_macc_ps(dz20,fscal,fjz0);
2288
2289             }
2290
2291             /**************************
2292              * CALCULATE INTERACTIONS *
2293              **************************/
2294
2295             if (gmx_mm_any_lt(rsq21,rcutoff2))
2296             {
2297
2298             r21              = _mm_mul_ps(rsq21,rinv21);
2299             r21              = _mm_andnot_ps(dummy_mask,r21);
2300
2301             /* EWALD ELECTROSTATICS */
2302
2303             /* Analytical PME correction */
2304             zeta2            = _mm_mul_ps(beta2,rsq21);
2305             rinv3            = _mm_mul_ps(rinvsq21,rinv21);
2306             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
2307             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
2308             felec            = _mm_mul_ps(qq21,felec);
2309
2310             cutoff_mask      = _mm_cmplt_ps(rsq21,rcutoff2);
2311
2312             fscal            = felec;
2313
2314             fscal            = _mm_and_ps(fscal,cutoff_mask);
2315
2316             fscal            = _mm_andnot_ps(dummy_mask,fscal);
2317
2318              /* Update vectorial force */
2319             fix2             = _mm_macc_ps(dx21,fscal,fix2);
2320             fiy2             = _mm_macc_ps(dy21,fscal,fiy2);
2321             fiz2             = _mm_macc_ps(dz21,fscal,fiz2);
2322
2323             fjx1             = _mm_macc_ps(dx21,fscal,fjx1);
2324             fjy1             = _mm_macc_ps(dy21,fscal,fjy1);
2325             fjz1             = _mm_macc_ps(dz21,fscal,fjz1);
2326
2327             }
2328
2329             /**************************
2330              * CALCULATE INTERACTIONS *
2331              **************************/
2332
2333             if (gmx_mm_any_lt(rsq22,rcutoff2))
2334             {
2335
2336             r22              = _mm_mul_ps(rsq22,rinv22);
2337             r22              = _mm_andnot_ps(dummy_mask,r22);
2338
2339             /* EWALD ELECTROSTATICS */
2340
2341             /* Analytical PME correction */
2342             zeta2            = _mm_mul_ps(beta2,rsq22);
2343             rinv3            = _mm_mul_ps(rinvsq22,rinv22);
2344             pmecorrF         = avx128fma_pmecorrF_f(zeta2);
2345             felec            = _mm_macc_ps(pmecorrF,beta3,rinv3);
2346             felec            = _mm_mul_ps(qq22,felec);
2347
2348             cutoff_mask      = _mm_cmplt_ps(rsq22,rcutoff2);
2349
2350             fscal            = felec;
2351
2352             fscal            = _mm_and_ps(fscal,cutoff_mask);
2353
2354             fscal            = _mm_andnot_ps(dummy_mask,fscal);
2355
2356              /* Update vectorial force */
2357             fix2             = _mm_macc_ps(dx22,fscal,fix2);
2358             fiy2             = _mm_macc_ps(dy22,fscal,fiy2);
2359             fiz2             = _mm_macc_ps(dz22,fscal,fiz2);
2360
2361             fjx2             = _mm_macc_ps(dx22,fscal,fjx2);
2362             fjy2             = _mm_macc_ps(dy22,fscal,fjy2);
2363             fjz2             = _mm_macc_ps(dz22,fscal,fjz2);
2364
2365             }
2366
2367             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
2368             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
2369             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
2370             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
2371
2372             gmx_mm_decrement_3rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,
2373                                                    fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,fjx2,fjy2,fjz2);
2374
2375             /* Inner loop uses 309 flops */
2376         }
2377
2378         /* End of innermost loop */
2379
2380         gmx_mm_update_iforce_3atom_swizzle_ps(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,
2381                                               f+i_coord_offset,fshift+i_shift_offset);
2382
2383         /* Increment number of inner iterations */
2384         inneriter                  += j_index_end - j_index_start;
2385
2386         /* Outer loop uses 18 flops */
2387     }
2388
2389     /* Increment number of outer iterations */
2390     outeriter        += nri;
2391
2392     /* Update outer/inner flops */
2393
2394     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W3W3_F,outeriter*18 + inneriter*309);
2395 }