src/gromacs/gmxlib/nonbonded/nb_kernel_sse2_single/nb_kernel_ElecEwSh_VdwLJSh_GeomW4W4_sse2_single.c

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