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

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