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

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