src/gromacs/gmxlib/nonbonded/nb_kernel_avx_128_fma_single/nb_kernel_ElecEwSh_VdwLJEwSh_GeomW4W4_avx_128_fma_single.c

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