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

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