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

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