src/gmxlib/nonbonded/nb_kernel_sse2_single/nb_kernel_ElecCSTab_VdwNone_GeomW3W3_sse2_single.c

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