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