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