src/gromacs/gmxlib/nonbonded/nb_kernel_sse4_1_single/nb_kernel_ElecEwSh_VdwLJSh_GeomW4W4_sse4_1_single.c

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