src/gromacs/gmxlib/nonbonded/nb_kernel_sse2_single/nb_kernel_ElecCSTab_VdwCSTab_GeomW3W3_sse2_single.c

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