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

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