src/gmxlib/nonbonded/nb_kernel_sse2_single/nb_kernel_ElecCSTab_VdwNone_GeomW4P1_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_GeomW4P1_VF_sse2_single
  38  * Electrostatics interaction: CubicSplineTable
  39  * VdW interaction:            None
  40  * Geometry:                   Water4-Particle
  41  * Calculate force/pot:        PotentialAndForce
  42  */
  43 void
  44 nb_kernel_ElecCSTab_VdwNone_GeomW4P1_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              j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
  62     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
  63     real             shX,shY,shZ,rcutoff_scalar;
  64     real             *shiftvec,*fshift,*x,*f;
  65     __m128           tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
  66     int              vdwioffset1;
  67     __m128           ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
  68     int              vdwioffset2;
  69     __m128           ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
  70     int              vdwioffset3;
  71     __m128           ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
  72     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
  73     __m128           jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
  74     __m128           dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
  75     __m128           dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
  76     __m128           dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
  77     __m128           velec,felec,velecsum,facel,crf,krf,krf2;
  78     real             *charge;
  79     __m128i          vfitab;
  80     __m128i          ifour       = _mm_set1_epi32(4);
  81     __m128           rt,vfeps,vftabscale,Y,F,G,H,Heps,Fp,VV,FF;
  82     real             *vftab;
  83     __m128           dummy_mask,cutoff_mask;
  84     __m128           signbit = _mm_castsi128_ps( _mm_set1_epi32(0x80000000) );
  85     __m128           one     = _mm_set1_ps(1.0);
  86     __m128           two     = _mm_set1_ps(2.0);
  87     x                = xx[0];
  88     f                = ff[0];
  89
  90     nri              = nlist->nri;
  91     iinr             = nlist->iinr;
  92     jindex           = nlist->jindex;
  93     jjnr             = nlist->jjnr;
  94     shiftidx         = nlist->shift;
  95     gid              = nlist->gid;
  96     shiftvec         = fr->shift_vec[0];
  97     fshift           = fr->fshift[0];
  98     facel            = _mm_set1_ps(fr->epsfac);
  99     charge           = mdatoms->chargeA;
 100
 101     vftab            = kernel_data->table_elec->data;
 102     vftabscale       = _mm_set1_ps(kernel_data->table_elec->scale);
 103
 104     /* Setup water-specific parameters */
 105     inr              = nlist->iinr[0];
 106     iq1              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+1]));
 107     iq2              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+2]));
 108     iq3              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+3]));
 109
 110     /* Avoid stupid compiler warnings */
 111     jnrA = jnrB = jnrC = jnrD = 0;
 112     j_coord_offsetA = 0;
 113     j_coord_offsetB = 0;
 114     j_coord_offsetC = 0;
 115     j_coord_offsetD = 0;
 116
 117     outeriter        = 0;
 118     inneriter        = 0;
 119
 120     /* Start outer loop over neighborlists */
 121     for(iidx=0; iidx<nri; iidx++)
 122     {
 123         /* Load shift vector for this list */
 124         i_shift_offset   = DIM*shiftidx[iidx];
 125         shX              = shiftvec[i_shift_offset+XX];
 126         shY              = shiftvec[i_shift_offset+YY];
 127         shZ              = shiftvec[i_shift_offset+ZZ];
 128
 129         /* Load limits for loop over neighbors */
 130         j_index_start    = jindex[iidx];
 131         j_index_end      = jindex[iidx+1];
 132
 133         /* Get outer coordinate index */
 134         inr              = iinr[iidx];
 135         i_coord_offset   = DIM*inr;
 136
 137         /* Load i particle coords and add shift vector */
 138         ix1              = _mm_set1_ps(shX + x[i_coord_offset+DIM*1+XX]);
 139         iy1              = _mm_set1_ps(shY + x[i_coord_offset+DIM*1+YY]);
 140         iz1              = _mm_set1_ps(shZ + x[i_coord_offset+DIM*1+ZZ]);
 141         ix2              = _mm_set1_ps(shX + x[i_coord_offset+DIM*2+XX]);
 142         iy2              = _mm_set1_ps(shY + x[i_coord_offset+DIM*2+YY]);
 143         iz2              = _mm_set1_ps(shZ + x[i_coord_offset+DIM*2+ZZ]);
 144         ix3              = _mm_set1_ps(shX + x[i_coord_offset+DIM*3+XX]);
 145         iy3              = _mm_set1_ps(shY + x[i_coord_offset+DIM*3+YY]);
 146         iz3              = _mm_set1_ps(shZ + x[i_coord_offset+DIM*3+ZZ]);
 147
 148         fix1             = _mm_setzero_ps();
 149         fiy1             = _mm_setzero_ps();
 150         fiz1             = _mm_setzero_ps();
 151         fix2             = _mm_setzero_ps();
 152         fiy2             = _mm_setzero_ps();
 153         fiz2             = _mm_setzero_ps();
 154         fix3             = _mm_setzero_ps();
 155         fiy3             = _mm_setzero_ps();
 156         fiz3             = _mm_setzero_ps();
 157
 158         /* Reset potential sums */
 159         velecsum         = _mm_setzero_ps();
 160
 161         /* Start inner kernel loop */
 162         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
 163         {
 164
 165             /* Get j neighbor index, and coordinate index */
 166             jnrA             = jjnr[jidx];
 167             jnrB             = jjnr[jidx+1];
 168             jnrC             = jjnr[jidx+2];
 169             jnrD             = jjnr[jidx+3];
 170
 171             j_coord_offsetA  = DIM*jnrA;
 172             j_coord_offsetB  = DIM*jnrB;
 173             j_coord_offsetC  = DIM*jnrC;
 174             j_coord_offsetD  = DIM*jnrD;
 175
 176             /* load j atom coordinates */
 177             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
 178                                               x+j_coord_offsetC,x+j_coord_offsetD,
 179                                               &jx0,&jy0,&jz0);
 180
 181             /* Calculate displacement vector */
 182             dx10             = _mm_sub_ps(ix1,jx0);
 183             dy10             = _mm_sub_ps(iy1,jy0);
 184             dz10             = _mm_sub_ps(iz1,jz0);
 185             dx20             = _mm_sub_ps(ix2,jx0);
 186             dy20             = _mm_sub_ps(iy2,jy0);
 187             dz20             = _mm_sub_ps(iz2,jz0);
 188             dx30             = _mm_sub_ps(ix3,jx0);
 189             dy30             = _mm_sub_ps(iy3,jy0);
 190             dz30             = _mm_sub_ps(iz3,jz0);
 191
 192             /* Calculate squared distance and things based on it */
 193             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
 194             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
 195             rsq30            = gmx_mm_calc_rsq_ps(dx30,dy30,dz30);
 196
 197             rinv10           = gmx_mm_invsqrt_ps(rsq10);
 198             rinv20           = gmx_mm_invsqrt_ps(rsq20);
 199             rinv30           = gmx_mm_invsqrt_ps(rsq30);
 200
 201             /* Load parameters for j particles */
 202             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
 203                                                               charge+jnrC+0,charge+jnrD+0);
 204
 205             /**************************
 206              * CALCULATE INTERACTIONS *
 207              **************************/
 208
 209             r10              = _mm_mul_ps(rsq10,rinv10);
 210
 211             /* Compute parameters for interactions between i and j atoms */
 212             qq10             = _mm_mul_ps(iq1,jq0);
 213
 214             /* Calculate table index by multiplying r with table scale and truncate to integer */
 215             rt               = _mm_mul_ps(r10,vftabscale);
 216             vfitab           = _mm_cvttps_epi32(rt);
 217             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
 218             vfitab           = _mm_slli_epi32(vfitab,2);
 219
 220             /* CUBIC SPLINE TABLE ELECTROSTATICS */
 221             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
 222             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
 223             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
 224             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
 225             _MM_TRANSPOSE4_PS(Y,F,G,H);
 226             Heps             = _mm_mul_ps(vfeps,H);
 227             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
 228             VV               = _mm_add_ps(Y,_mm_mul_ps(vfeps,Fp));
 229             velec            = _mm_mul_ps(qq10,VV);
 230             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
 231             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq10,FF),_mm_mul_ps(vftabscale,rinv10)));
 232
 233             /* Update potential sum for this i atom from the interaction with this j atom. */
 234             velecsum         = _mm_add_ps(velecsum,velec);
 235
 236             fscal            = felec;
 237
 238             /* Calculate temporary vectorial force */
 239             tx               = _mm_mul_ps(fscal,dx10);
 240             ty               = _mm_mul_ps(fscal,dy10);
 241             tz               = _mm_mul_ps(fscal,dz10);
 242
 243             /* Update vectorial force */
 244             fix1             = _mm_add_ps(fix1,tx);
 245             fiy1             = _mm_add_ps(fiy1,ty);
 246             fiz1             = _mm_add_ps(fiz1,tz);
 247
 248             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(f+j_coord_offsetA,f+j_coord_offsetB,
 249                                                    f+j_coord_offsetC,f+j_coord_offsetD,
 250                                                    tx,ty,tz);
 251
 252             /**************************
 253              * CALCULATE INTERACTIONS *
 254              **************************/
 255
 256             r20              = _mm_mul_ps(rsq20,rinv20);
 257
 258             /* Compute parameters for interactions between i and j atoms */
 259             qq20             = _mm_mul_ps(iq2,jq0);
 260
 261             /* Calculate table index by multiplying r with table scale and truncate to integer */
 262             rt               = _mm_mul_ps(r20,vftabscale);
 263             vfitab           = _mm_cvttps_epi32(rt);
 264             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
 265             vfitab           = _mm_slli_epi32(vfitab,2);
 266
 267             /* CUBIC SPLINE TABLE ELECTROSTATICS */
 268             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
 269             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
 270             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
 271             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
 272             _MM_TRANSPOSE4_PS(Y,F,G,H);
 273             Heps             = _mm_mul_ps(vfeps,H);
 274             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
 275             VV               = _mm_add_ps(Y,_mm_mul_ps(vfeps,Fp));
 276             velec            = _mm_mul_ps(qq20,VV);
 277             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
 278             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq20,FF),_mm_mul_ps(vftabscale,rinv20)));
 279
 280             /* Update potential sum for this i atom from the interaction with this j atom. */
 281             velecsum         = _mm_add_ps(velecsum,velec);
 282
 283             fscal            = felec;
 284
 285             /* Calculate temporary vectorial force */
 286             tx               = _mm_mul_ps(fscal,dx20);
 287             ty               = _mm_mul_ps(fscal,dy20);
 288             tz               = _mm_mul_ps(fscal,dz20);
 289
 290             /* Update vectorial force */
 291             fix2             = _mm_add_ps(fix2,tx);
 292             fiy2             = _mm_add_ps(fiy2,ty);
 293             fiz2             = _mm_add_ps(fiz2,tz);
 294
 295             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(f+j_coord_offsetA,f+j_coord_offsetB,
 296                                                    f+j_coord_offsetC,f+j_coord_offsetD,
 297                                                    tx,ty,tz);
 298
 299             /**************************
 300              * CALCULATE INTERACTIONS *
 301              **************************/
 302
 303             r30              = _mm_mul_ps(rsq30,rinv30);
 304
 305             /* Compute parameters for interactions between i and j atoms */
 306             qq30             = _mm_mul_ps(iq3,jq0);
 307
 308             /* Calculate table index by multiplying r with table scale and truncate to integer */
 309             rt               = _mm_mul_ps(r30,vftabscale);
 310             vfitab           = _mm_cvttps_epi32(rt);
 311             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
 312             vfitab           = _mm_slli_epi32(vfitab,2);
 313
 314             /* CUBIC SPLINE TABLE ELECTROSTATICS */
 315             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
 316             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
 317             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
 318             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
 319             _MM_TRANSPOSE4_PS(Y,F,G,H);
 320             Heps             = _mm_mul_ps(vfeps,H);
 321             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
 322             VV               = _mm_add_ps(Y,_mm_mul_ps(vfeps,Fp));
 323             velec            = _mm_mul_ps(qq30,VV);
 324             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
 325             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq30,FF),_mm_mul_ps(vftabscale,rinv30)));
 326
 327             /* Update potential sum for this i atom from the interaction with this j atom. */
 328             velecsum         = _mm_add_ps(velecsum,velec);
 329
 330             fscal            = felec;
 331
 332             /* Calculate temporary vectorial force */
 333             tx               = _mm_mul_ps(fscal,dx30);
 334             ty               = _mm_mul_ps(fscal,dy30);
 335             tz               = _mm_mul_ps(fscal,dz30);
 336
 337             /* Update vectorial force */
 338             fix3             = _mm_add_ps(fix3,tx);
 339             fiy3             = _mm_add_ps(fiy3,ty);
 340             fiz3             = _mm_add_ps(fiz3,tz);
 341
 342             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(f+j_coord_offsetA,f+j_coord_offsetB,
 343                                                    f+j_coord_offsetC,f+j_coord_offsetD,
 344                                                    tx,ty,tz);
 345
 346             /* Inner loop uses 129 flops */
 347         }
 348
 349         if(jidx<j_index_end)
 350         {
 351
 352             /* Get j neighbor index, and coordinate index */
 353             jnrA             = jjnr[jidx];
 354             jnrB             = jjnr[jidx+1];
 355             jnrC             = jjnr[jidx+2];
 356             jnrD             = jjnr[jidx+3];
 357
 358             /* Sign of each element will be negative for non-real atoms.
 359              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
 360              * so use it as val = _mm_andnot_ps(mask,val) to clear dummy entries.
 361              */
 362             dummy_mask = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
 363             jnrA       = (jnrA>=0) ? jnrA : 0;
 364             jnrB       = (jnrB>=0) ? jnrB : 0;
 365             jnrC       = (jnrC>=0) ? jnrC : 0;
 366             jnrD       = (jnrD>=0) ? jnrD : 0;
 367
 368             j_coord_offsetA  = DIM*jnrA;
 369             j_coord_offsetB  = DIM*jnrB;
 370             j_coord_offsetC  = DIM*jnrC;
 371             j_coord_offsetD  = DIM*jnrD;
 372
 373             /* load j atom coordinates */
 374             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
 375                                               x+j_coord_offsetC,x+j_coord_offsetD,
 376                                               &jx0,&jy0,&jz0);
 377
 378             /* Calculate displacement vector */
 379             dx10             = _mm_sub_ps(ix1,jx0);
 380             dy10             = _mm_sub_ps(iy1,jy0);
 381             dz10             = _mm_sub_ps(iz1,jz0);
 382             dx20             = _mm_sub_ps(ix2,jx0);
 383             dy20             = _mm_sub_ps(iy2,jy0);
 384             dz20             = _mm_sub_ps(iz2,jz0);
 385             dx30             = _mm_sub_ps(ix3,jx0);
 386             dy30             = _mm_sub_ps(iy3,jy0);
 387             dz30             = _mm_sub_ps(iz3,jz0);
 388
 389             /* Calculate squared distance and things based on it */
 390             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
 391             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
 392             rsq30            = gmx_mm_calc_rsq_ps(dx30,dy30,dz30);
 393
 394             rinv10           = gmx_mm_invsqrt_ps(rsq10);
 395             rinv20           = gmx_mm_invsqrt_ps(rsq20);
 396             rinv30           = gmx_mm_invsqrt_ps(rsq30);
 397
 398             /* Load parameters for j particles */
 399             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
 400                                                               charge+jnrC+0,charge+jnrD+0);
 401
 402             /**************************
 403              * CALCULATE INTERACTIONS *
 404              **************************/
 405
 406             r10              = _mm_mul_ps(rsq10,rinv10);
 407             r10              = _mm_andnot_ps(dummy_mask,r10);
 408
 409             /* Compute parameters for interactions between i and j atoms */
 410             qq10             = _mm_mul_ps(iq1,jq0);
 411
 412             /* Calculate table index by multiplying r with table scale and truncate to integer */
 413             rt               = _mm_mul_ps(r10,vftabscale);
 414             vfitab           = _mm_cvttps_epi32(rt);
 415             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
 416             vfitab           = _mm_slli_epi32(vfitab,2);
 417
 418             /* CUBIC SPLINE TABLE ELECTROSTATICS */
 419             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
 420             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
 421             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
 422             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
 423             _MM_TRANSPOSE4_PS(Y,F,G,H);
 424             Heps             = _mm_mul_ps(vfeps,H);
 425             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
 426             VV               = _mm_add_ps(Y,_mm_mul_ps(vfeps,Fp));
 427             velec            = _mm_mul_ps(qq10,VV);
 428             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
 429             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq10,FF),_mm_mul_ps(vftabscale,rinv10)));
 430
 431             /* Update potential sum for this i atom from the interaction with this j atom. */
 432             velec            = _mm_andnot_ps(dummy_mask,velec);
 433             velecsum         = _mm_add_ps(velecsum,velec);
 434
 435             fscal            = felec;
 436
 437             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 438
 439             /* Calculate temporary vectorial force */
 440             tx               = _mm_mul_ps(fscal,dx10);
 441             ty               = _mm_mul_ps(fscal,dy10);
 442             tz               = _mm_mul_ps(fscal,dz10);
 443
 444             /* Update vectorial force */
 445             fix1             = _mm_add_ps(fix1,tx);
 446             fiy1             = _mm_add_ps(fiy1,ty);
 447             fiz1             = _mm_add_ps(fiz1,tz);
 448
 449             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(f+j_coord_offsetA,f+j_coord_offsetB,
 450                                                    f+j_coord_offsetC,f+j_coord_offsetD,
 451                                                    tx,ty,tz);
 452
 453             /**************************
 454              * CALCULATE INTERACTIONS *
 455              **************************/
 456
 457             r20              = _mm_mul_ps(rsq20,rinv20);
 458             r20              = _mm_andnot_ps(dummy_mask,r20);
 459
 460             /* Compute parameters for interactions between i and j atoms */
 461             qq20             = _mm_mul_ps(iq2,jq0);
 462
 463             /* Calculate table index by multiplying r with table scale and truncate to integer */
 464             rt               = _mm_mul_ps(r20,vftabscale);
 465             vfitab           = _mm_cvttps_epi32(rt);
 466             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
 467             vfitab           = _mm_slli_epi32(vfitab,2);
 468
 469             /* CUBIC SPLINE TABLE ELECTROSTATICS */
 470             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
 471             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
 472             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
 473             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
 474             _MM_TRANSPOSE4_PS(Y,F,G,H);
 475             Heps             = _mm_mul_ps(vfeps,H);
 476             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
 477             VV               = _mm_add_ps(Y,_mm_mul_ps(vfeps,Fp));
 478             velec            = _mm_mul_ps(qq20,VV);
 479             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
 480             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq20,FF),_mm_mul_ps(vftabscale,rinv20)));
 481
 482             /* Update potential sum for this i atom from the interaction with this j atom. */
 483             velec            = _mm_andnot_ps(dummy_mask,velec);
 484             velecsum         = _mm_add_ps(velecsum,velec);
 485
 486             fscal            = felec;
 487
 488             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 489
 490             /* Calculate temporary vectorial force */
 491             tx               = _mm_mul_ps(fscal,dx20);
 492             ty               = _mm_mul_ps(fscal,dy20);
 493             tz               = _mm_mul_ps(fscal,dz20);
 494
 495             /* Update vectorial force */
 496             fix2             = _mm_add_ps(fix2,tx);
 497             fiy2             = _mm_add_ps(fiy2,ty);
 498             fiz2             = _mm_add_ps(fiz2,tz);
 499
 500             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(f+j_coord_offsetA,f+j_coord_offsetB,
 501                                                    f+j_coord_offsetC,f+j_coord_offsetD,
 502                                                    tx,ty,tz);
 503
 504             /**************************
 505              * CALCULATE INTERACTIONS *
 506              **************************/
 507
 508             r30              = _mm_mul_ps(rsq30,rinv30);
 509             r30              = _mm_andnot_ps(dummy_mask,r30);
 510
 511             /* Compute parameters for interactions between i and j atoms */
 512             qq30             = _mm_mul_ps(iq3,jq0);
 513
 514             /* Calculate table index by multiplying r with table scale and truncate to integer */
 515             rt               = _mm_mul_ps(r30,vftabscale);
 516             vfitab           = _mm_cvttps_epi32(rt);
 517             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
 518             vfitab           = _mm_slli_epi32(vfitab,2);
 519
 520             /* CUBIC SPLINE TABLE ELECTROSTATICS */
 521             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
 522             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
 523             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
 524             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
 525             _MM_TRANSPOSE4_PS(Y,F,G,H);
 526             Heps             = _mm_mul_ps(vfeps,H);
 527             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
 528             VV               = _mm_add_ps(Y,_mm_mul_ps(vfeps,Fp));
 529             velec            = _mm_mul_ps(qq30,VV);
 530             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
 531             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq30,FF),_mm_mul_ps(vftabscale,rinv30)));
 532
 533             /* Update potential sum for this i atom from the interaction with this j atom. */
 534             velec            = _mm_andnot_ps(dummy_mask,velec);
 535             velecsum         = _mm_add_ps(velecsum,velec);
 536
 537             fscal            = felec;
 538
 539             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 540
 541             /* Calculate temporary vectorial force */
 542             tx               = _mm_mul_ps(fscal,dx30);
 543             ty               = _mm_mul_ps(fscal,dy30);
 544             tz               = _mm_mul_ps(fscal,dz30);
 545
 546             /* Update vectorial force */
 547             fix3             = _mm_add_ps(fix3,tx);
 548             fiy3             = _mm_add_ps(fiy3,ty);
 549             fiz3             = _mm_add_ps(fiz3,tz);
 550
 551             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(f+j_coord_offsetA,f+j_coord_offsetB,
 552                                                    f+j_coord_offsetC,f+j_coord_offsetD,
 553                                                    tx,ty,tz);
 554
 555             /* Inner loop uses 132 flops */
 556         }
 557
 558         /* End of innermost loop */
 559
 560         gmx_mm_update_iforce_3atom_swizzle_ps(fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
 561                                               f+i_coord_offset+DIM,fshift+i_shift_offset);
 562
 563         ggid                        = gid[iidx];
 564         /* Update potential energies */
 565         gmx_mm_update_1pot_ps(velecsum,kernel_data->energygrp_elec+ggid);
 566
 567         /* Increment number of inner iterations */
 568         inneriter                  += j_index_end - j_index_start;
 569
 570         /* Outer loop uses 28 flops */
 571     }
 572
 573     /* Increment number of outer iterations */
 574     outeriter        += nri;
 575
 576     /* Update outer/inner flops */
 577
 578     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_W4_VF,outeriter*28 + inneriter*132);
 579 }
 580 /*
 581  * Gromacs nonbonded kernel:   nb_kernel_ElecCSTab_VdwNone_GeomW4P1_F_sse2_single
 582  * Electrostatics interaction: CubicSplineTable
 583  * VdW interaction:            None
 584  * Geometry:                   Water4-Particle
 585  * Calculate force/pot:        Force
 586  */
 587 void
 588 nb_kernel_ElecCSTab_VdwNone_GeomW4P1_F_sse2_single
 589                     (t_nblist * gmx_restrict                nlist,
 590                      rvec * gmx_restrict                    xx,
 591                      rvec * gmx_restrict                    ff,
 592                      t_forcerec * gmx_restrict              fr,
 593                      t_mdatoms * gmx_restrict               mdatoms,
 594                      nb_kernel_data_t * gmx_restrict        kernel_data,
 595                      t_nrnb * gmx_restrict                  nrnb)
 596 {
 597     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
 598      * just 0 for non-waters.
 599      * Suffixes A,B,C,D refer to j loop unrolling done with SSE, e.g. for the four different
 600      * jnr indices corresponding to data put in the four positions in the SIMD register.
 601      */
 602     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
 603     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
 604     int              jnrA,jnrB,jnrC,jnrD;
 605     int              j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
 606     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
 607     real             shX,shY,shZ,rcutoff_scalar;
 608     real             *shiftvec,*fshift,*x,*f;
 609     __m128           tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
 610     int              vdwioffset1;
 611     __m128           ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
 612     int              vdwioffset2;
 613     __m128           ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
 614     int              vdwioffset3;
 615     __m128           ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
 616     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
 617     __m128           jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
 618     __m128           dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
 619     __m128           dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
 620     __m128           dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
 621     __m128           velec,felec,velecsum,facel,crf,krf,krf2;
 622     real             *charge;
 623     __m128i          vfitab;
 624     __m128i          ifour       = _mm_set1_epi32(4);
 625     __m128           rt,vfeps,vftabscale,Y,F,G,H,Heps,Fp,VV,FF;
 626     real             *vftab;
 627     __m128           dummy_mask,cutoff_mask;
 628     __m128           signbit = _mm_castsi128_ps( _mm_set1_epi32(0x80000000) );
 629     __m128           one     = _mm_set1_ps(1.0);
 630     __m128           two     = _mm_set1_ps(2.0);
 631     x                = xx[0];
 632     f                = ff[0];
 633
 634     nri              = nlist->nri;
 635     iinr             = nlist->iinr;
 636     jindex           = nlist->jindex;
 637     jjnr             = nlist->jjnr;
 638     shiftidx         = nlist->shift;
 639     gid              = nlist->gid;
 640     shiftvec         = fr->shift_vec[0];
 641     fshift           = fr->fshift[0];
 642     facel            = _mm_set1_ps(fr->epsfac);
 643     charge           = mdatoms->chargeA;
 644
 645     vftab            = kernel_data->table_elec->data;
 646     vftabscale       = _mm_set1_ps(kernel_data->table_elec->scale);
 647
 648     /* Setup water-specific parameters */
 649     inr              = nlist->iinr[0];
 650     iq1              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+1]));
 651     iq2              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+2]));
 652     iq3              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+3]));
 653
 654     /* Avoid stupid compiler warnings */
 655     jnrA = jnrB = jnrC = jnrD = 0;
 656     j_coord_offsetA = 0;
 657     j_coord_offsetB = 0;
 658     j_coord_offsetC = 0;
 659     j_coord_offsetD = 0;
 660
 661     outeriter        = 0;
 662     inneriter        = 0;
 663
 664     /* Start outer loop over neighborlists */
 665     for(iidx=0; iidx<nri; iidx++)
 666     {
 667         /* Load shift vector for this list */
 668         i_shift_offset   = DIM*shiftidx[iidx];
 669         shX              = shiftvec[i_shift_offset+XX];
 670         shY              = shiftvec[i_shift_offset+YY];
 671         shZ              = shiftvec[i_shift_offset+ZZ];
 672
 673         /* Load limits for loop over neighbors */
 674         j_index_start    = jindex[iidx];
 675         j_index_end      = jindex[iidx+1];
 676
 677         /* Get outer coordinate index */
 678         inr              = iinr[iidx];
 679         i_coord_offset   = DIM*inr;
 680
 681         /* Load i particle coords and add shift vector */
 682         ix1              = _mm_set1_ps(shX + x[i_coord_offset+DIM*1+XX]);
 683         iy1              = _mm_set1_ps(shY + x[i_coord_offset+DIM*1+YY]);
 684         iz1              = _mm_set1_ps(shZ + x[i_coord_offset+DIM*1+ZZ]);
 685         ix2              = _mm_set1_ps(shX + x[i_coord_offset+DIM*2+XX]);
 686         iy2              = _mm_set1_ps(shY + x[i_coord_offset+DIM*2+YY]);
 687         iz2              = _mm_set1_ps(shZ + x[i_coord_offset+DIM*2+ZZ]);
 688         ix3              = _mm_set1_ps(shX + x[i_coord_offset+DIM*3+XX]);
 689         iy3              = _mm_set1_ps(shY + x[i_coord_offset+DIM*3+YY]);
 690         iz3              = _mm_set1_ps(shZ + x[i_coord_offset+DIM*3+ZZ]);
 691
 692         fix1             = _mm_setzero_ps();
 693         fiy1             = _mm_setzero_ps();
 694         fiz1             = _mm_setzero_ps();
 695         fix2             = _mm_setzero_ps();
 696         fiy2             = _mm_setzero_ps();
 697         fiz2             = _mm_setzero_ps();
 698         fix3             = _mm_setzero_ps();
 699         fiy3             = _mm_setzero_ps();
 700         fiz3             = _mm_setzero_ps();
 701
 702         /* Start inner kernel loop */
 703         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
 704         {
 705
 706             /* Get j neighbor index, and coordinate index */
 707             jnrA             = jjnr[jidx];
 708             jnrB             = jjnr[jidx+1];
 709             jnrC             = jjnr[jidx+2];
 710             jnrD             = jjnr[jidx+3];
 711
 712             j_coord_offsetA  = DIM*jnrA;
 713             j_coord_offsetB  = DIM*jnrB;
 714             j_coord_offsetC  = DIM*jnrC;
 715             j_coord_offsetD  = DIM*jnrD;
 716
 717             /* load j atom coordinates */
 718             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
 719                                               x+j_coord_offsetC,x+j_coord_offsetD,
 720                                               &jx0,&jy0,&jz0);
 721
 722             /* Calculate displacement vector */
 723             dx10             = _mm_sub_ps(ix1,jx0);
 724             dy10             = _mm_sub_ps(iy1,jy0);
 725             dz10             = _mm_sub_ps(iz1,jz0);
 726             dx20             = _mm_sub_ps(ix2,jx0);
 727             dy20             = _mm_sub_ps(iy2,jy0);
 728             dz20             = _mm_sub_ps(iz2,jz0);
 729             dx30             = _mm_sub_ps(ix3,jx0);
 730             dy30             = _mm_sub_ps(iy3,jy0);
 731             dz30             = _mm_sub_ps(iz3,jz0);
 732
 733             /* Calculate squared distance and things based on it */
 734             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
 735             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
 736             rsq30            = gmx_mm_calc_rsq_ps(dx30,dy30,dz30);
 737
 738             rinv10           = gmx_mm_invsqrt_ps(rsq10);
 739             rinv20           = gmx_mm_invsqrt_ps(rsq20);
 740             rinv30           = gmx_mm_invsqrt_ps(rsq30);
 741
 742             /* Load parameters for j particles */
 743             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
 744                                                               charge+jnrC+0,charge+jnrD+0);
 745
 746             /**************************
 747              * CALCULATE INTERACTIONS *
 748              **************************/
 749
 750             r10              = _mm_mul_ps(rsq10,rinv10);
 751
 752             /* Compute parameters for interactions between i and j atoms */
 753             qq10             = _mm_mul_ps(iq1,jq0);
 754
 755             /* Calculate table index by multiplying r with table scale and truncate to integer */
 756             rt               = _mm_mul_ps(r10,vftabscale);
 757             vfitab           = _mm_cvttps_epi32(rt);
 758             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
 759             vfitab           = _mm_slli_epi32(vfitab,2);
 760
 761             /* CUBIC SPLINE TABLE ELECTROSTATICS */
 762             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
 763             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
 764             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
 765             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
 766             _MM_TRANSPOSE4_PS(Y,F,G,H);
 767             Heps             = _mm_mul_ps(vfeps,H);
 768             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
 769             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
 770             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq10,FF),_mm_mul_ps(vftabscale,rinv10)));
 771
 772             fscal            = felec;
 773
 774             /* Calculate temporary vectorial force */
 775             tx               = _mm_mul_ps(fscal,dx10);
 776             ty               = _mm_mul_ps(fscal,dy10);
 777             tz               = _mm_mul_ps(fscal,dz10);
 778
 779             /* Update vectorial force */
 780             fix1             = _mm_add_ps(fix1,tx);
 781             fiy1             = _mm_add_ps(fiy1,ty);
 782             fiz1             = _mm_add_ps(fiz1,tz);
 783
 784             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(f+j_coord_offsetA,f+j_coord_offsetB,
 785                                                    f+j_coord_offsetC,f+j_coord_offsetD,
 786                                                    tx,ty,tz);
 787
 788             /**************************
 789              * CALCULATE INTERACTIONS *
 790              **************************/
 791
 792             r20              = _mm_mul_ps(rsq20,rinv20);
 793
 794             /* Compute parameters for interactions between i and j atoms */
 795             qq20             = _mm_mul_ps(iq2,jq0);
 796
 797             /* Calculate table index by multiplying r with table scale and truncate to integer */
 798             rt               = _mm_mul_ps(r20,vftabscale);
 799             vfitab           = _mm_cvttps_epi32(rt);
 800             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
 801             vfitab           = _mm_slli_epi32(vfitab,2);
 802
 803             /* CUBIC SPLINE TABLE ELECTROSTATICS */
 804             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
 805             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
 806             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
 807             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
 808             _MM_TRANSPOSE4_PS(Y,F,G,H);
 809             Heps             = _mm_mul_ps(vfeps,H);
 810             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
 811             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
 812             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq20,FF),_mm_mul_ps(vftabscale,rinv20)));
 813
 814             fscal            = felec;
 815
 816             /* Calculate temporary vectorial force */
 817             tx               = _mm_mul_ps(fscal,dx20);
 818             ty               = _mm_mul_ps(fscal,dy20);
 819             tz               = _mm_mul_ps(fscal,dz20);
 820
 821             /* Update vectorial force */
 822             fix2             = _mm_add_ps(fix2,tx);
 823             fiy2             = _mm_add_ps(fiy2,ty);
 824             fiz2             = _mm_add_ps(fiz2,tz);
 825
 826             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(f+j_coord_offsetA,f+j_coord_offsetB,
 827                                                    f+j_coord_offsetC,f+j_coord_offsetD,
 828                                                    tx,ty,tz);
 829
 830             /**************************
 831              * CALCULATE INTERACTIONS *
 832              **************************/
 833
 834             r30              = _mm_mul_ps(rsq30,rinv30);
 835
 836             /* Compute parameters for interactions between i and j atoms */
 837             qq30             = _mm_mul_ps(iq3,jq0);
 838
 839             /* Calculate table index by multiplying r with table scale and truncate to integer */
 840             rt               = _mm_mul_ps(r30,vftabscale);
 841             vfitab           = _mm_cvttps_epi32(rt);
 842             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
 843             vfitab           = _mm_slli_epi32(vfitab,2);
 844
 845             /* CUBIC SPLINE TABLE ELECTROSTATICS */
 846             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
 847             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
 848             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
 849             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
 850             _MM_TRANSPOSE4_PS(Y,F,G,H);
 851             Heps             = _mm_mul_ps(vfeps,H);
 852             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
 853             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
 854             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq30,FF),_mm_mul_ps(vftabscale,rinv30)));
 855
 856             fscal            = felec;
 857
 858             /* Calculate temporary vectorial force */
 859             tx               = _mm_mul_ps(fscal,dx30);
 860             ty               = _mm_mul_ps(fscal,dy30);
 861             tz               = _mm_mul_ps(fscal,dz30);
 862
 863             /* Update vectorial force */
 864             fix3             = _mm_add_ps(fix3,tx);
 865             fiy3             = _mm_add_ps(fiy3,ty);
 866             fiz3             = _mm_add_ps(fiz3,tz);
 867
 868             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(f+j_coord_offsetA,f+j_coord_offsetB,
 869                                                    f+j_coord_offsetC,f+j_coord_offsetD,
 870                                                    tx,ty,tz);
 871
 872             /* Inner loop uses 117 flops */
 873         }
 874
 875         if(jidx<j_index_end)
 876         {
 877
 878             /* Get j neighbor index, and coordinate index */
 879             jnrA             = jjnr[jidx];
 880             jnrB             = jjnr[jidx+1];
 881             jnrC             = jjnr[jidx+2];
 882             jnrD             = jjnr[jidx+3];
 883
 884             /* Sign of each element will be negative for non-real atoms.
 885              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
 886              * so use it as val = _mm_andnot_ps(mask,val) to clear dummy entries.
 887              */
 888             dummy_mask = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
 889             jnrA       = (jnrA>=0) ? jnrA : 0;
 890             jnrB       = (jnrB>=0) ? jnrB : 0;
 891             jnrC       = (jnrC>=0) ? jnrC : 0;
 892             jnrD       = (jnrD>=0) ? jnrD : 0;
 893
 894             j_coord_offsetA  = DIM*jnrA;
 895             j_coord_offsetB  = DIM*jnrB;
 896             j_coord_offsetC  = DIM*jnrC;
 897             j_coord_offsetD  = DIM*jnrD;
 898
 899             /* load j atom coordinates */
 900             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
 901                                               x+j_coord_offsetC,x+j_coord_offsetD,
 902                                               &jx0,&jy0,&jz0);
 903
 904             /* Calculate displacement vector */
 905             dx10             = _mm_sub_ps(ix1,jx0);
 906             dy10             = _mm_sub_ps(iy1,jy0);
 907             dz10             = _mm_sub_ps(iz1,jz0);
 908             dx20             = _mm_sub_ps(ix2,jx0);
 909             dy20             = _mm_sub_ps(iy2,jy0);
 910             dz20             = _mm_sub_ps(iz2,jz0);
 911             dx30             = _mm_sub_ps(ix3,jx0);
 912             dy30             = _mm_sub_ps(iy3,jy0);
 913             dz30             = _mm_sub_ps(iz3,jz0);
 914
 915             /* Calculate squared distance and things based on it */
 916             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
 917             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
 918             rsq30            = gmx_mm_calc_rsq_ps(dx30,dy30,dz30);
 919
 920             rinv10           = gmx_mm_invsqrt_ps(rsq10);
 921             rinv20           = gmx_mm_invsqrt_ps(rsq20);
 922             rinv30           = gmx_mm_invsqrt_ps(rsq30);
 923
 924             /* Load parameters for j particles */
 925             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
 926                                                               charge+jnrC+0,charge+jnrD+0);
 927
 928             /**************************
 929              * CALCULATE INTERACTIONS *
 930              **************************/
 931
 932             r10              = _mm_mul_ps(rsq10,rinv10);
 933             r10              = _mm_andnot_ps(dummy_mask,r10);
 934
 935             /* Compute parameters for interactions between i and j atoms */
 936             qq10             = _mm_mul_ps(iq1,jq0);
 937
 938             /* Calculate table index by multiplying r with table scale and truncate to integer */
 939             rt               = _mm_mul_ps(r10,vftabscale);
 940             vfitab           = _mm_cvttps_epi32(rt);
 941             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
 942             vfitab           = _mm_slli_epi32(vfitab,2);
 943
 944             /* CUBIC SPLINE TABLE ELECTROSTATICS */
 945             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
 946             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
 947             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
 948             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
 949             _MM_TRANSPOSE4_PS(Y,F,G,H);
 950             Heps             = _mm_mul_ps(vfeps,H);
 951             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
 952             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
 953             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq10,FF),_mm_mul_ps(vftabscale,rinv10)));
 954
 955             fscal            = felec;
 956
 957             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 958
 959             /* Calculate temporary vectorial force */
 960             tx               = _mm_mul_ps(fscal,dx10);
 961             ty               = _mm_mul_ps(fscal,dy10);
 962             tz               = _mm_mul_ps(fscal,dz10);
 963
 964             /* Update vectorial force */
 965             fix1             = _mm_add_ps(fix1,tx);
 966             fiy1             = _mm_add_ps(fiy1,ty);
 967             fiz1             = _mm_add_ps(fiz1,tz);
 968
 969             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(f+j_coord_offsetA,f+j_coord_offsetB,
 970                                                    f+j_coord_offsetC,f+j_coord_offsetD,
 971                                                    tx,ty,tz);
 972
 973             /**************************
 974              * CALCULATE INTERACTIONS *
 975              **************************/
 976
 977             r20              = _mm_mul_ps(rsq20,rinv20);
 978             r20              = _mm_andnot_ps(dummy_mask,r20);
 979
 980             /* Compute parameters for interactions between i and j atoms */
 981             qq20             = _mm_mul_ps(iq2,jq0);
 982
 983             /* Calculate table index by multiplying r with table scale and truncate to integer */
 984             rt               = _mm_mul_ps(r20,vftabscale);
 985             vfitab           = _mm_cvttps_epi32(rt);
 986             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
 987             vfitab           = _mm_slli_epi32(vfitab,2);
 988
 989             /* CUBIC SPLINE TABLE ELECTROSTATICS */
 990             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
 991             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
 992             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
 993             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
 994             _MM_TRANSPOSE4_PS(Y,F,G,H);
 995             Heps             = _mm_mul_ps(vfeps,H);
 996             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
 997             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
 998             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq20,FF),_mm_mul_ps(vftabscale,rinv20)));
 999
1000             fscal            = felec;
1001
1002             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1003
1004             /* Calculate temporary vectorial force */
1005             tx               = _mm_mul_ps(fscal,dx20);
1006             ty               = _mm_mul_ps(fscal,dy20);
1007             tz               = _mm_mul_ps(fscal,dz20);
1008
1009             /* Update vectorial force */
1010             fix2             = _mm_add_ps(fix2,tx);
1011             fiy2             = _mm_add_ps(fiy2,ty);
1012             fiz2             = _mm_add_ps(fiz2,tz);
1013
1014             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(f+j_coord_offsetA,f+j_coord_offsetB,
1015                                                    f+j_coord_offsetC,f+j_coord_offsetD,
1016                                                    tx,ty,tz);
1017
1018             /**************************
1019              * CALCULATE INTERACTIONS *
1020              **************************/
1021
1022             r30              = _mm_mul_ps(rsq30,rinv30);
1023             r30              = _mm_andnot_ps(dummy_mask,r30);
1024
1025             /* Compute parameters for interactions between i and j atoms */
1026             qq30             = _mm_mul_ps(iq3,jq0);
1027
1028             /* Calculate table index by multiplying r with table scale and truncate to integer */
1029             rt               = _mm_mul_ps(r30,vftabscale);
1030             vfitab           = _mm_cvttps_epi32(rt);
1031             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
1032             vfitab           = _mm_slli_epi32(vfitab,2);
1033
1034             /* CUBIC SPLINE TABLE ELECTROSTATICS */
1035             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
1036             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
1037             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
1038             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
1039             _MM_TRANSPOSE4_PS(Y,F,G,H);
1040             Heps             = _mm_mul_ps(vfeps,H);
1041             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
1042             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
1043             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq30,FF),_mm_mul_ps(vftabscale,rinv30)));
1044
1045             fscal            = felec;
1046
1047             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1048
1049             /* Calculate temporary vectorial force */
1050             tx               = _mm_mul_ps(fscal,dx30);
1051             ty               = _mm_mul_ps(fscal,dy30);
1052             tz               = _mm_mul_ps(fscal,dz30);
1053
1054             /* Update vectorial force */
1055             fix3             = _mm_add_ps(fix3,tx);
1056             fiy3             = _mm_add_ps(fiy3,ty);
1057             fiz3             = _mm_add_ps(fiz3,tz);
1058
1059             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(f+j_coord_offsetA,f+j_coord_offsetB,
1060                                                    f+j_coord_offsetC,f+j_coord_offsetD,
1061                                                    tx,ty,tz);
1062
1063             /* Inner loop uses 120 flops */
1064         }
1065
1066         /* End of innermost loop */
1067
1068         gmx_mm_update_iforce_3atom_swizzle_ps(fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
1069                                               f+i_coord_offset+DIM,fshift+i_shift_offset);
1070
1071         /* Increment number of inner iterations */
1072         inneriter                  += j_index_end - j_index_start;
1073
1074         /* Outer loop uses 27 flops */
1075     }
1076
1077     /* Increment number of outer iterations */
1078     outeriter        += nri;
1079
1080     /* Update outer/inner flops */
1081
1082     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_W4_F,outeriter*27 + inneriter*120);
1083 }