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