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