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