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