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

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