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