src/gmxlib/nonbonded/nb_kernel_avx_128_fma_double/nb_kernel_ElecEw_VdwCSTab_GeomW4W4_avx_128_fma_double.c

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