src/gmxlib/nonbonded/nb_kernel_sse4_1_single/nb_kernel_ElecRF_VdwCSTab_GeomW4W4_sse4_1_single.c

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