src/gmxlib/nonbonded/nb_kernel_sse4_1_single/nb_kernel_ElecEw_VdwNone_GeomW4P1_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_ElecEw_VdwNone_GeomW4P1_VF_sse4_1_single
  38  * Electrostatics interaction: Ewald
  39  * VdW interaction:            None
  40  * Geometry:                   Water4-Particle
  41  * Calculate force/pot:        PotentialAndForce
  42  */
  43 void
  44 nb_kernel_ElecEw_VdwNone_GeomW4P1_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              vdwioffset1;
  70     __m128           ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
  71     int              vdwioffset2;
  72     __m128           ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
  73     int              vdwioffset3;
  74     __m128           ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
  75     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
  76     __m128           jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
  77     __m128           dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
  78     __m128           dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
  79     __m128           dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
  80     __m128           velec,felec,velecsum,facel,crf,krf,krf2;
  81     real             *charge;
  82     __m128i          ewitab;
  83     __m128           ewtabscale,eweps,sh_ewald,ewrt,ewtabhalfspace,ewtabF,ewtabFn,ewtabD,ewtabV;
  84     real             *ewtab;
  85     __m128           dummy_mask,cutoff_mask;
  86     __m128           signbit = _mm_castsi128_ps( _mm_set1_epi32(0x80000000) );
  87     __m128           one     = _mm_set1_ps(1.0);
  88     __m128           two     = _mm_set1_ps(2.0);
  89     x                = xx[0];
  90     f                = ff[0];
  91
  92     nri              = nlist->nri;
  93     iinr             = nlist->iinr;
  94     jindex           = nlist->jindex;
  95     jjnr             = nlist->jjnr;
  96     shiftidx         = nlist->shift;
  97     gid              = nlist->gid;
  98     shiftvec         = fr->shift_vec[0];
  99     fshift           = fr->fshift[0];
 100     facel            = _mm_set1_ps(fr->epsfac);
 101     charge           = mdatoms->chargeA;
 102
 103     sh_ewald         = _mm_set1_ps(fr->ic->sh_ewald);
 104     ewtab            = fr->ic->tabq_coul_FDV0;
 105     ewtabscale       = _mm_set1_ps(fr->ic->tabq_scale);
 106     ewtabhalfspace   = _mm_set1_ps(0.5/fr->ic->tabq_scale);
 107
 108     /* Setup water-specific parameters */
 109     inr              = nlist->iinr[0];
 110     iq1              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+1]));
 111     iq2              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+2]));
 112     iq3              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+3]));
 113
 114     /* Avoid stupid compiler warnings */
 115     jnrA = jnrB = jnrC = jnrD = 0;
 116     j_coord_offsetA = 0;
 117     j_coord_offsetB = 0;
 118     j_coord_offsetC = 0;
 119     j_coord_offsetD = 0;
 120
 121     outeriter        = 0;
 122     inneriter        = 0;
 123
 124     for(iidx=0;iidx<4*DIM;iidx++)
 125     {
 126         scratch[iidx] = 0.0;
 127     }
 128
 129     /* Start outer loop over neighborlists */
 130     for(iidx=0; iidx<nri; iidx++)
 131     {
 132         /* Load shift vector for this list */
 133         i_shift_offset   = DIM*shiftidx[iidx];
 134
 135         /* Load limits for loop over neighbors */
 136         j_index_start    = jindex[iidx];
 137         j_index_end      = jindex[iidx+1];
 138
 139         /* Get outer coordinate index */
 140         inr              = iinr[iidx];
 141         i_coord_offset   = DIM*inr;
 142
 143         /* Load i particle coords and add shift vector */
 144         gmx_mm_load_shift_and_3rvec_broadcast_ps(shiftvec+i_shift_offset,x+i_coord_offset+DIM,
 145                                                  &ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
 146
 147         fix1             = _mm_setzero_ps();
 148         fiy1             = _mm_setzero_ps();
 149         fiz1             = _mm_setzero_ps();
 150         fix2             = _mm_setzero_ps();
 151         fiy2             = _mm_setzero_ps();
 152         fiz2             = _mm_setzero_ps();
 153         fix3             = _mm_setzero_ps();
 154         fiy3             = _mm_setzero_ps();
 155         fiz3             = _mm_setzero_ps();
 156
 157         /* Reset potential sums */
 158         velecsum         = _mm_setzero_ps();
 159
 160         /* Start inner kernel loop */
 161         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
 162         {
 163
 164             /* Get j neighbor index, and coordinate index */
 165             jnrA             = jjnr[jidx];
 166             jnrB             = jjnr[jidx+1];
 167             jnrC             = jjnr[jidx+2];
 168             jnrD             = jjnr[jidx+3];
 169             j_coord_offsetA  = DIM*jnrA;
 170             j_coord_offsetB  = DIM*jnrB;
 171             j_coord_offsetC  = DIM*jnrC;
 172             j_coord_offsetD  = DIM*jnrD;
 173
 174             /* load j atom coordinates */
 175             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
 176                                               x+j_coord_offsetC,x+j_coord_offsetD,
 177                                               &jx0,&jy0,&jz0);
 178
 179             /* Calculate displacement vector */
 180             dx10             = _mm_sub_ps(ix1,jx0);
 181             dy10             = _mm_sub_ps(iy1,jy0);
 182             dz10             = _mm_sub_ps(iz1,jz0);
 183             dx20             = _mm_sub_ps(ix2,jx0);
 184             dy20             = _mm_sub_ps(iy2,jy0);
 185             dz20             = _mm_sub_ps(iz2,jz0);
 186             dx30             = _mm_sub_ps(ix3,jx0);
 187             dy30             = _mm_sub_ps(iy3,jy0);
 188             dz30             = _mm_sub_ps(iz3,jz0);
 189
 190             /* Calculate squared distance and things based on it */
 191             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
 192             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
 193             rsq30            = gmx_mm_calc_rsq_ps(dx30,dy30,dz30);
 194
 195             rinv10           = gmx_mm_invsqrt_ps(rsq10);
 196             rinv20           = gmx_mm_invsqrt_ps(rsq20);
 197             rinv30           = gmx_mm_invsqrt_ps(rsq30);
 198
 199             rinvsq10         = _mm_mul_ps(rinv10,rinv10);
 200             rinvsq20         = _mm_mul_ps(rinv20,rinv20);
 201             rinvsq30         = _mm_mul_ps(rinv30,rinv30);
 202
 203             /* Load parameters for j particles */
 204             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
 205                                                               charge+jnrC+0,charge+jnrD+0);
 206
 207             /**************************
 208              * CALCULATE INTERACTIONS *
 209              **************************/
 210
 211             r10              = _mm_mul_ps(rsq10,rinv10);
 212
 213             /* Compute parameters for interactions between i and j atoms */
 214             qq10             = _mm_mul_ps(iq1,jq0);
 215
 216             /* EWALD ELECTROSTATICS */
 217
 218             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 219             ewrt             = _mm_mul_ps(r10,ewtabscale);
 220             ewitab           = _mm_cvttps_epi32(ewrt);
 221             eweps            = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
 222             ewitab           = _mm_slli_epi32(ewitab,2);
 223             ewtabF           = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
 224             ewtabD           = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
 225             ewtabV           = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
 226             ewtabFn          = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
 227             _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
 228             felec            = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
 229             velec            = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
 230             velec            = _mm_mul_ps(qq10,_mm_sub_ps(rinv10,velec));
 231             felec            = _mm_mul_ps(_mm_mul_ps(qq10,rinv10),_mm_sub_ps(rinvsq10,felec));
 232
 233             /* Update potential sum for this i atom from the interaction with this j atom. */
 234             velecsum         = _mm_add_ps(velecsum,velec);
 235
 236             fscal            = felec;
 237
 238             /* Calculate temporary vectorial force */
 239             tx               = _mm_mul_ps(fscal,dx10);
 240             ty               = _mm_mul_ps(fscal,dy10);
 241             tz               = _mm_mul_ps(fscal,dz10);
 242
 243             /* Update vectorial force */
 244             fix1             = _mm_add_ps(fix1,tx);
 245             fiy1             = _mm_add_ps(fiy1,ty);
 246             fiz1             = _mm_add_ps(fiz1,tz);
 247
 248             fjptrA             = f+j_coord_offsetA;
 249             fjptrB             = f+j_coord_offsetB;
 250             fjptrC             = f+j_coord_offsetC;
 251             fjptrD             = f+j_coord_offsetD;
 252             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
 253
 254             /**************************
 255              * CALCULATE INTERACTIONS *
 256              **************************/
 257
 258             r20              = _mm_mul_ps(rsq20,rinv20);
 259
 260             /* Compute parameters for interactions between i and j atoms */
 261             qq20             = _mm_mul_ps(iq2,jq0);
 262
 263             /* EWALD ELECTROSTATICS */
 264
 265             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 266             ewrt             = _mm_mul_ps(r20,ewtabscale);
 267             ewitab           = _mm_cvttps_epi32(ewrt);
 268             eweps            = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
 269             ewitab           = _mm_slli_epi32(ewitab,2);
 270             ewtabF           = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
 271             ewtabD           = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
 272             ewtabV           = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
 273             ewtabFn          = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
 274             _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
 275             felec            = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
 276             velec            = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
 277             velec            = _mm_mul_ps(qq20,_mm_sub_ps(rinv20,velec));
 278             felec            = _mm_mul_ps(_mm_mul_ps(qq20,rinv20),_mm_sub_ps(rinvsq20,felec));
 279
 280             /* Update potential sum for this i atom from the interaction with this j atom. */
 281             velecsum         = _mm_add_ps(velecsum,velec);
 282
 283             fscal            = felec;
 284
 285             /* Calculate temporary vectorial force */
 286             tx               = _mm_mul_ps(fscal,dx20);
 287             ty               = _mm_mul_ps(fscal,dy20);
 288             tz               = _mm_mul_ps(fscal,dz20);
 289
 290             /* Update vectorial force */
 291             fix2             = _mm_add_ps(fix2,tx);
 292             fiy2             = _mm_add_ps(fiy2,ty);
 293             fiz2             = _mm_add_ps(fiz2,tz);
 294
 295             fjptrA             = f+j_coord_offsetA;
 296             fjptrB             = f+j_coord_offsetB;
 297             fjptrC             = f+j_coord_offsetC;
 298             fjptrD             = f+j_coord_offsetD;
 299             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
 300
 301             /**************************
 302              * CALCULATE INTERACTIONS *
 303              **************************/
 304
 305             r30              = _mm_mul_ps(rsq30,rinv30);
 306
 307             /* Compute parameters for interactions between i and j atoms */
 308             qq30             = _mm_mul_ps(iq3,jq0);
 309
 310             /* EWALD ELECTROSTATICS */
 311
 312             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 313             ewrt             = _mm_mul_ps(r30,ewtabscale);
 314             ewitab           = _mm_cvttps_epi32(ewrt);
 315             eweps            = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
 316             ewitab           = _mm_slli_epi32(ewitab,2);
 317             ewtabF           = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
 318             ewtabD           = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
 319             ewtabV           = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
 320             ewtabFn          = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
 321             _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
 322             felec            = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
 323             velec            = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
 324             velec            = _mm_mul_ps(qq30,_mm_sub_ps(rinv30,velec));
 325             felec            = _mm_mul_ps(_mm_mul_ps(qq30,rinv30),_mm_sub_ps(rinvsq30,felec));
 326
 327             /* Update potential sum for this i atom from the interaction with this j atom. */
 328             velecsum         = _mm_add_ps(velecsum,velec);
 329
 330             fscal            = felec;
 331
 332             /* Calculate temporary vectorial force */
 333             tx               = _mm_mul_ps(fscal,dx30);
 334             ty               = _mm_mul_ps(fscal,dy30);
 335             tz               = _mm_mul_ps(fscal,dz30);
 336
 337             /* Update vectorial force */
 338             fix3             = _mm_add_ps(fix3,tx);
 339             fiy3             = _mm_add_ps(fiy3,ty);
 340             fiz3             = _mm_add_ps(fiz3,tz);
 341
 342             fjptrA             = f+j_coord_offsetA;
 343             fjptrB             = f+j_coord_offsetB;
 344             fjptrC             = f+j_coord_offsetC;
 345             fjptrD             = f+j_coord_offsetD;
 346             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
 347
 348             /* Inner loop uses 123 flops */
 349         }
 350
 351         if(jidx<j_index_end)
 352         {
 353
 354             /* Get j neighbor index, and coordinate index */
 355             jnrlistA         = jjnr[jidx];
 356             jnrlistB         = jjnr[jidx+1];
 357             jnrlistC         = jjnr[jidx+2];
 358             jnrlistD         = jjnr[jidx+3];
 359             /* Sign of each element will be negative for non-real atoms.
 360              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
 361              * so use it as val = _mm_andnot_ps(mask,val) to clear dummy entries.
 362              */
 363             dummy_mask = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
 364             jnrA       = (jnrlistA>=0) ? jnrlistA : 0;
 365             jnrB       = (jnrlistB>=0) ? jnrlistB : 0;
 366             jnrC       = (jnrlistC>=0) ? jnrlistC : 0;
 367             jnrD       = (jnrlistD>=0) ? jnrlistD : 0;
 368             j_coord_offsetA  = DIM*jnrA;
 369             j_coord_offsetB  = DIM*jnrB;
 370             j_coord_offsetC  = DIM*jnrC;
 371             j_coord_offsetD  = DIM*jnrD;
 372
 373             /* load j atom coordinates */
 374             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
 375                                               x+j_coord_offsetC,x+j_coord_offsetD,
 376                                               &jx0,&jy0,&jz0);
 377
 378             /* Calculate displacement vector */
 379             dx10             = _mm_sub_ps(ix1,jx0);
 380             dy10             = _mm_sub_ps(iy1,jy0);
 381             dz10             = _mm_sub_ps(iz1,jz0);
 382             dx20             = _mm_sub_ps(ix2,jx0);
 383             dy20             = _mm_sub_ps(iy2,jy0);
 384             dz20             = _mm_sub_ps(iz2,jz0);
 385             dx30             = _mm_sub_ps(ix3,jx0);
 386             dy30             = _mm_sub_ps(iy3,jy0);
 387             dz30             = _mm_sub_ps(iz3,jz0);
 388
 389             /* Calculate squared distance and things based on it */
 390             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
 391             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
 392             rsq30            = gmx_mm_calc_rsq_ps(dx30,dy30,dz30);
 393
 394             rinv10           = gmx_mm_invsqrt_ps(rsq10);
 395             rinv20           = gmx_mm_invsqrt_ps(rsq20);
 396             rinv30           = gmx_mm_invsqrt_ps(rsq30);
 397
 398             rinvsq10         = _mm_mul_ps(rinv10,rinv10);
 399             rinvsq20         = _mm_mul_ps(rinv20,rinv20);
 400             rinvsq30         = _mm_mul_ps(rinv30,rinv30);
 401
 402             /* Load parameters for j particles */
 403             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
 404                                                               charge+jnrC+0,charge+jnrD+0);
 405
 406             /**************************
 407              * CALCULATE INTERACTIONS *
 408              **************************/
 409
 410             r10              = _mm_mul_ps(rsq10,rinv10);
 411             r10              = _mm_andnot_ps(dummy_mask,r10);
 412
 413             /* Compute parameters for interactions between i and j atoms */
 414             qq10             = _mm_mul_ps(iq1,jq0);
 415
 416             /* EWALD ELECTROSTATICS */
 417
 418             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 419             ewrt             = _mm_mul_ps(r10,ewtabscale);
 420             ewitab           = _mm_cvttps_epi32(ewrt);
 421             eweps            = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
 422             ewitab           = _mm_slli_epi32(ewitab,2);
 423             ewtabF           = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
 424             ewtabD           = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
 425             ewtabV           = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
 426             ewtabFn          = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
 427             _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
 428             felec            = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
 429             velec            = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
 430             velec            = _mm_mul_ps(qq10,_mm_sub_ps(rinv10,velec));
 431             felec            = _mm_mul_ps(_mm_mul_ps(qq10,rinv10),_mm_sub_ps(rinvsq10,felec));
 432
 433             /* Update potential sum for this i atom from the interaction with this j atom. */
 434             velec            = _mm_andnot_ps(dummy_mask,velec);
 435             velecsum         = _mm_add_ps(velecsum,velec);
 436
 437             fscal            = felec;
 438
 439             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 440
 441             /* Calculate temporary vectorial force */
 442             tx               = _mm_mul_ps(fscal,dx10);
 443             ty               = _mm_mul_ps(fscal,dy10);
 444             tz               = _mm_mul_ps(fscal,dz10);
 445
 446             /* Update vectorial force */
 447             fix1             = _mm_add_ps(fix1,tx);
 448             fiy1             = _mm_add_ps(fiy1,ty);
 449             fiz1             = _mm_add_ps(fiz1,tz);
 450
 451             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
 452             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
 453             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
 454             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
 455             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
 456
 457             /**************************
 458              * CALCULATE INTERACTIONS *
 459              **************************/
 460
 461             r20              = _mm_mul_ps(rsq20,rinv20);
 462             r20              = _mm_andnot_ps(dummy_mask,r20);
 463
 464             /* Compute parameters for interactions between i and j atoms */
 465             qq20             = _mm_mul_ps(iq2,jq0);
 466
 467             /* EWALD ELECTROSTATICS */
 468
 469             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 470             ewrt             = _mm_mul_ps(r20,ewtabscale);
 471             ewitab           = _mm_cvttps_epi32(ewrt);
 472             eweps            = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
 473             ewitab           = _mm_slli_epi32(ewitab,2);
 474             ewtabF           = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
 475             ewtabD           = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
 476             ewtabV           = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
 477             ewtabFn          = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
 478             _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
 479             felec            = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
 480             velec            = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
 481             velec            = _mm_mul_ps(qq20,_mm_sub_ps(rinv20,velec));
 482             felec            = _mm_mul_ps(_mm_mul_ps(qq20,rinv20),_mm_sub_ps(rinvsq20,felec));
 483
 484             /* Update potential sum for this i atom from the interaction with this j atom. */
 485             velec            = _mm_andnot_ps(dummy_mask,velec);
 486             velecsum         = _mm_add_ps(velecsum,velec);
 487
 488             fscal            = felec;
 489
 490             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 491
 492             /* Calculate temporary vectorial force */
 493             tx               = _mm_mul_ps(fscal,dx20);
 494             ty               = _mm_mul_ps(fscal,dy20);
 495             tz               = _mm_mul_ps(fscal,dz20);
 496
 497             /* Update vectorial force */
 498             fix2             = _mm_add_ps(fix2,tx);
 499             fiy2             = _mm_add_ps(fiy2,ty);
 500             fiz2             = _mm_add_ps(fiz2,tz);
 501
 502             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
 503             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
 504             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
 505             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
 506             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
 507
 508             /**************************
 509              * CALCULATE INTERACTIONS *
 510              **************************/
 511
 512             r30              = _mm_mul_ps(rsq30,rinv30);
 513             r30              = _mm_andnot_ps(dummy_mask,r30);
 514
 515             /* Compute parameters for interactions between i and j atoms */
 516             qq30             = _mm_mul_ps(iq3,jq0);
 517
 518             /* EWALD ELECTROSTATICS */
 519
 520             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 521             ewrt             = _mm_mul_ps(r30,ewtabscale);
 522             ewitab           = _mm_cvttps_epi32(ewrt);
 523             eweps            = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
 524             ewitab           = _mm_slli_epi32(ewitab,2);
 525             ewtabF           = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
 526             ewtabD           = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
 527             ewtabV           = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
 528             ewtabFn          = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
 529             _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
 530             felec            = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
 531             velec            = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
 532             velec            = _mm_mul_ps(qq30,_mm_sub_ps(rinv30,velec));
 533             felec            = _mm_mul_ps(_mm_mul_ps(qq30,rinv30),_mm_sub_ps(rinvsq30,felec));
 534
 535             /* Update potential sum for this i atom from the interaction with this j atom. */
 536             velec            = _mm_andnot_ps(dummy_mask,velec);
 537             velecsum         = _mm_add_ps(velecsum,velec);
 538
 539             fscal            = felec;
 540
 541             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 542
 543             /* Calculate temporary vectorial force */
 544             tx               = _mm_mul_ps(fscal,dx30);
 545             ty               = _mm_mul_ps(fscal,dy30);
 546             tz               = _mm_mul_ps(fscal,dz30);
 547
 548             /* Update vectorial force */
 549             fix3             = _mm_add_ps(fix3,tx);
 550             fiy3             = _mm_add_ps(fiy3,ty);
 551             fiz3             = _mm_add_ps(fiz3,tz);
 552
 553             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
 554             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
 555             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
 556             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
 557             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
 558
 559             /* Inner loop uses 126 flops */
 560         }
 561
 562         /* End of innermost loop */
 563
 564         gmx_mm_update_iforce_3atom_swizzle_ps(fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
 565                                               f+i_coord_offset+DIM,fshift+i_shift_offset);
 566
 567         ggid                        = gid[iidx];
 568         /* Update potential energies */
 569         gmx_mm_update_1pot_ps(velecsum,kernel_data->energygrp_elec+ggid);
 570
 571         /* Increment number of inner iterations */
 572         inneriter                  += j_index_end - j_index_start;
 573
 574         /* Outer loop uses 19 flops */
 575     }
 576
 577     /* Increment number of outer iterations */
 578     outeriter        += nri;
 579
 580     /* Update outer/inner flops */
 581
 582     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_W4_VF,outeriter*19 + inneriter*126);
 583 }
 584 /*
 585  * Gromacs nonbonded kernel:   nb_kernel_ElecEw_VdwNone_GeomW4P1_F_sse4_1_single
 586  * Electrostatics interaction: Ewald
 587  * VdW interaction:            None
 588  * Geometry:                   Water4-Particle
 589  * Calculate force/pot:        Force
 590  */
 591 void
 592 nb_kernel_ElecEw_VdwNone_GeomW4P1_F_sse4_1_single
 593                     (t_nblist * gmx_restrict                nlist,
 594                      rvec * gmx_restrict                    xx,
 595                      rvec * gmx_restrict                    ff,
 596                      t_forcerec * gmx_restrict              fr,
 597                      t_mdatoms * gmx_restrict               mdatoms,
 598                      nb_kernel_data_t * gmx_restrict        kernel_data,
 599                      t_nrnb * gmx_restrict                  nrnb)
 600 {
 601     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
 602      * just 0 for non-waters.
 603      * Suffixes A,B,C,D refer to j loop unrolling done with SSE, e.g. for the four different
 604      * jnr indices corresponding to data put in the four positions in the SIMD register.
 605      */
 606     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
 607     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
 608     int              jnrA,jnrB,jnrC,jnrD;
 609     int              jnrlistA,jnrlistB,jnrlistC,jnrlistD;
 610     int              j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
 611     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
 612     real             rcutoff_scalar;
 613     real             *shiftvec,*fshift,*x,*f;
 614     real             *fjptrA,*fjptrB,*fjptrC,*fjptrD;
 615     real             scratch[4*DIM];
 616     __m128           tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
 617     int              vdwioffset1;
 618     __m128           ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
 619     int              vdwioffset2;
 620     __m128           ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
 621     int              vdwioffset3;
 622     __m128           ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
 623     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
 624     __m128           jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
 625     __m128           dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
 626     __m128           dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
 627     __m128           dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
 628     __m128           velec,felec,velecsum,facel,crf,krf,krf2;
 629     real             *charge;
 630     __m128i          ewitab;
 631     __m128           ewtabscale,eweps,sh_ewald,ewrt,ewtabhalfspace,ewtabF,ewtabFn,ewtabD,ewtabV;
 632     real             *ewtab;
 633     __m128           dummy_mask,cutoff_mask;
 634     __m128           signbit = _mm_castsi128_ps( _mm_set1_epi32(0x80000000) );
 635     __m128           one     = _mm_set1_ps(1.0);
 636     __m128           two     = _mm_set1_ps(2.0);
 637     x                = xx[0];
 638     f                = ff[0];
 639
 640     nri              = nlist->nri;
 641     iinr             = nlist->iinr;
 642     jindex           = nlist->jindex;
 643     jjnr             = nlist->jjnr;
 644     shiftidx         = nlist->shift;
 645     gid              = nlist->gid;
 646     shiftvec         = fr->shift_vec[0];
 647     fshift           = fr->fshift[0];
 648     facel            = _mm_set1_ps(fr->epsfac);
 649     charge           = mdatoms->chargeA;
 650
 651     sh_ewald         = _mm_set1_ps(fr->ic->sh_ewald);
 652     ewtab            = fr->ic->tabq_coul_F;
 653     ewtabscale       = _mm_set1_ps(fr->ic->tabq_scale);
 654     ewtabhalfspace   = _mm_set1_ps(0.5/fr->ic->tabq_scale);
 655
 656     /* Setup water-specific parameters */
 657     inr              = nlist->iinr[0];
 658     iq1              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+1]));
 659     iq2              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+2]));
 660     iq3              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+3]));
 661
 662     /* Avoid stupid compiler warnings */
 663     jnrA = jnrB = jnrC = jnrD = 0;
 664     j_coord_offsetA = 0;
 665     j_coord_offsetB = 0;
 666     j_coord_offsetC = 0;
 667     j_coord_offsetD = 0;
 668
 669     outeriter        = 0;
 670     inneriter        = 0;
 671
 672     for(iidx=0;iidx<4*DIM;iidx++)
 673     {
 674         scratch[iidx] = 0.0;
 675     }
 676
 677     /* Start outer loop over neighborlists */
 678     for(iidx=0; iidx<nri; iidx++)
 679     {
 680         /* Load shift vector for this list */
 681         i_shift_offset   = DIM*shiftidx[iidx];
 682
 683         /* Load limits for loop over neighbors */
 684         j_index_start    = jindex[iidx];
 685         j_index_end      = jindex[iidx+1];
 686
 687         /* Get outer coordinate index */
 688         inr              = iinr[iidx];
 689         i_coord_offset   = DIM*inr;
 690
 691         /* Load i particle coords and add shift vector */
 692         gmx_mm_load_shift_and_3rvec_broadcast_ps(shiftvec+i_shift_offset,x+i_coord_offset+DIM,
 693                                                  &ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
 694
 695         fix1             = _mm_setzero_ps();
 696         fiy1             = _mm_setzero_ps();
 697         fiz1             = _mm_setzero_ps();
 698         fix2             = _mm_setzero_ps();
 699         fiy2             = _mm_setzero_ps();
 700         fiz2             = _mm_setzero_ps();
 701         fix3             = _mm_setzero_ps();
 702         fiy3             = _mm_setzero_ps();
 703         fiz3             = _mm_setzero_ps();
 704
 705         /* Start inner kernel loop */
 706         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
 707         {
 708
 709             /* Get j neighbor index, and coordinate index */
 710             jnrA             = jjnr[jidx];
 711             jnrB             = jjnr[jidx+1];
 712             jnrC             = jjnr[jidx+2];
 713             jnrD             = jjnr[jidx+3];
 714             j_coord_offsetA  = DIM*jnrA;
 715             j_coord_offsetB  = DIM*jnrB;
 716             j_coord_offsetC  = DIM*jnrC;
 717             j_coord_offsetD  = DIM*jnrD;
 718
 719             /* load j atom coordinates */
 720             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
 721                                               x+j_coord_offsetC,x+j_coord_offsetD,
 722                                               &jx0,&jy0,&jz0);
 723
 724             /* Calculate displacement vector */
 725             dx10             = _mm_sub_ps(ix1,jx0);
 726             dy10             = _mm_sub_ps(iy1,jy0);
 727             dz10             = _mm_sub_ps(iz1,jz0);
 728             dx20             = _mm_sub_ps(ix2,jx0);
 729             dy20             = _mm_sub_ps(iy2,jy0);
 730             dz20             = _mm_sub_ps(iz2,jz0);
 731             dx30             = _mm_sub_ps(ix3,jx0);
 732             dy30             = _mm_sub_ps(iy3,jy0);
 733             dz30             = _mm_sub_ps(iz3,jz0);
 734
 735             /* Calculate squared distance and things based on it */
 736             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
 737             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
 738             rsq30            = gmx_mm_calc_rsq_ps(dx30,dy30,dz30);
 739
 740             rinv10           = gmx_mm_invsqrt_ps(rsq10);
 741             rinv20           = gmx_mm_invsqrt_ps(rsq20);
 742             rinv30           = gmx_mm_invsqrt_ps(rsq30);
 743
 744             rinvsq10         = _mm_mul_ps(rinv10,rinv10);
 745             rinvsq20         = _mm_mul_ps(rinv20,rinv20);
 746             rinvsq30         = _mm_mul_ps(rinv30,rinv30);
 747
 748             /* Load parameters for j particles */
 749             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
 750                                                               charge+jnrC+0,charge+jnrD+0);
 751
 752             /**************************
 753              * CALCULATE INTERACTIONS *
 754              **************************/
 755
 756             r10              = _mm_mul_ps(rsq10,rinv10);
 757
 758             /* Compute parameters for interactions between i and j atoms */
 759             qq10             = _mm_mul_ps(iq1,jq0);
 760
 761             /* EWALD ELECTROSTATICS */
 762
 763             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 764             ewrt             = _mm_mul_ps(r10,ewtabscale);
 765             ewitab           = _mm_cvttps_epi32(ewrt);
 766             eweps            = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
 767             gmx_mm_load_4pair_swizzle_ps(ewtab + gmx_mm_extract_epi32(ewitab,0),ewtab + gmx_mm_extract_epi32(ewitab,1),
 768                                          ewtab + gmx_mm_extract_epi32(ewitab,2),ewtab + gmx_mm_extract_epi32(ewitab,3),
 769                                          &ewtabF,&ewtabFn);
 770             felec            = _mm_add_ps(_mm_mul_ps( _mm_sub_ps(one,eweps),ewtabF),_mm_mul_ps(eweps,ewtabFn));
 771             felec            = _mm_mul_ps(_mm_mul_ps(qq10,rinv10),_mm_sub_ps(rinvsq10,felec));
 772
 773             fscal            = felec;
 774
 775             /* Calculate temporary vectorial force */
 776             tx               = _mm_mul_ps(fscal,dx10);
 777             ty               = _mm_mul_ps(fscal,dy10);
 778             tz               = _mm_mul_ps(fscal,dz10);
 779
 780             /* Update vectorial force */
 781             fix1             = _mm_add_ps(fix1,tx);
 782             fiy1             = _mm_add_ps(fiy1,ty);
 783             fiz1             = _mm_add_ps(fiz1,tz);
 784
 785             fjptrA             = f+j_coord_offsetA;
 786             fjptrB             = f+j_coord_offsetB;
 787             fjptrC             = f+j_coord_offsetC;
 788             fjptrD             = f+j_coord_offsetD;
 789             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
 790
 791             /**************************
 792              * CALCULATE INTERACTIONS *
 793              **************************/
 794
 795             r20              = _mm_mul_ps(rsq20,rinv20);
 796
 797             /* Compute parameters for interactions between i and j atoms */
 798             qq20             = _mm_mul_ps(iq2,jq0);
 799
 800             /* EWALD ELECTROSTATICS */
 801
 802             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 803             ewrt             = _mm_mul_ps(r20,ewtabscale);
 804             ewitab           = _mm_cvttps_epi32(ewrt);
 805             eweps            = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
 806             gmx_mm_load_4pair_swizzle_ps(ewtab + gmx_mm_extract_epi32(ewitab,0),ewtab + gmx_mm_extract_epi32(ewitab,1),
 807                                          ewtab + gmx_mm_extract_epi32(ewitab,2),ewtab + gmx_mm_extract_epi32(ewitab,3),
 808                                          &ewtabF,&ewtabFn);
 809             felec            = _mm_add_ps(_mm_mul_ps( _mm_sub_ps(one,eweps),ewtabF),_mm_mul_ps(eweps,ewtabFn));
 810             felec            = _mm_mul_ps(_mm_mul_ps(qq20,rinv20),_mm_sub_ps(rinvsq20,felec));
 811
 812             fscal            = felec;
 813
 814             /* Calculate temporary vectorial force */
 815             tx               = _mm_mul_ps(fscal,dx20);
 816             ty               = _mm_mul_ps(fscal,dy20);
 817             tz               = _mm_mul_ps(fscal,dz20);
 818
 819             /* Update vectorial force */
 820             fix2             = _mm_add_ps(fix2,tx);
 821             fiy2             = _mm_add_ps(fiy2,ty);
 822             fiz2             = _mm_add_ps(fiz2,tz);
 823
 824             fjptrA             = f+j_coord_offsetA;
 825             fjptrB             = f+j_coord_offsetB;
 826             fjptrC             = f+j_coord_offsetC;
 827             fjptrD             = f+j_coord_offsetD;
 828             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
 829
 830             /**************************
 831              * CALCULATE INTERACTIONS *
 832              **************************/
 833
 834             r30              = _mm_mul_ps(rsq30,rinv30);
 835
 836             /* Compute parameters for interactions between i and j atoms */
 837             qq30             = _mm_mul_ps(iq3,jq0);
 838
 839             /* EWALD ELECTROSTATICS */
 840
 841             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 842             ewrt             = _mm_mul_ps(r30,ewtabscale);
 843             ewitab           = _mm_cvttps_epi32(ewrt);
 844             eweps            = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
 845             gmx_mm_load_4pair_swizzle_ps(ewtab + gmx_mm_extract_epi32(ewitab,0),ewtab + gmx_mm_extract_epi32(ewitab,1),
 846                                          ewtab + gmx_mm_extract_epi32(ewitab,2),ewtab + gmx_mm_extract_epi32(ewitab,3),
 847                                          &ewtabF,&ewtabFn);
 848             felec            = _mm_add_ps(_mm_mul_ps( _mm_sub_ps(one,eweps),ewtabF),_mm_mul_ps(eweps,ewtabFn));
 849             felec            = _mm_mul_ps(_mm_mul_ps(qq30,rinv30),_mm_sub_ps(rinvsq30,felec));
 850
 851             fscal            = felec;
 852
 853             /* Calculate temporary vectorial force */
 854             tx               = _mm_mul_ps(fscal,dx30);
 855             ty               = _mm_mul_ps(fscal,dy30);
 856             tz               = _mm_mul_ps(fscal,dz30);
 857
 858             /* Update vectorial force */
 859             fix3             = _mm_add_ps(fix3,tx);
 860             fiy3             = _mm_add_ps(fiy3,ty);
 861             fiz3             = _mm_add_ps(fiz3,tz);
 862
 863             fjptrA             = f+j_coord_offsetA;
 864             fjptrB             = f+j_coord_offsetB;
 865             fjptrC             = f+j_coord_offsetC;
 866             fjptrD             = f+j_coord_offsetD;
 867             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
 868
 869             /* Inner loop uses 108 flops */
 870         }
 871
 872         if(jidx<j_index_end)
 873         {
 874
 875             /* Get j neighbor index, and coordinate index */
 876             jnrlistA         = jjnr[jidx];
 877             jnrlistB         = jjnr[jidx+1];
 878             jnrlistC         = jjnr[jidx+2];
 879             jnrlistD         = jjnr[jidx+3];
 880             /* Sign of each element will be negative for non-real atoms.
 881              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
 882              * so use it as val = _mm_andnot_ps(mask,val) to clear dummy entries.
 883              */
 884             dummy_mask = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
 885             jnrA       = (jnrlistA>=0) ? jnrlistA : 0;
 886             jnrB       = (jnrlistB>=0) ? jnrlistB : 0;
 887             jnrC       = (jnrlistC>=0) ? jnrlistC : 0;
 888             jnrD       = (jnrlistD>=0) ? jnrlistD : 0;
 889             j_coord_offsetA  = DIM*jnrA;
 890             j_coord_offsetB  = DIM*jnrB;
 891             j_coord_offsetC  = DIM*jnrC;
 892             j_coord_offsetD  = DIM*jnrD;
 893
 894             /* load j atom coordinates */
 895             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
 896                                               x+j_coord_offsetC,x+j_coord_offsetD,
 897                                               &jx0,&jy0,&jz0);
 898
 899             /* Calculate displacement vector */
 900             dx10             = _mm_sub_ps(ix1,jx0);
 901             dy10             = _mm_sub_ps(iy1,jy0);
 902             dz10             = _mm_sub_ps(iz1,jz0);
 903             dx20             = _mm_sub_ps(ix2,jx0);
 904             dy20             = _mm_sub_ps(iy2,jy0);
 905             dz20             = _mm_sub_ps(iz2,jz0);
 906             dx30             = _mm_sub_ps(ix3,jx0);
 907             dy30             = _mm_sub_ps(iy3,jy0);
 908             dz30             = _mm_sub_ps(iz3,jz0);
 909
 910             /* Calculate squared distance and things based on it */
 911             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
 912             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
 913             rsq30            = gmx_mm_calc_rsq_ps(dx30,dy30,dz30);
 914
 915             rinv10           = gmx_mm_invsqrt_ps(rsq10);
 916             rinv20           = gmx_mm_invsqrt_ps(rsq20);
 917             rinv30           = gmx_mm_invsqrt_ps(rsq30);
 918
 919             rinvsq10         = _mm_mul_ps(rinv10,rinv10);
 920             rinvsq20         = _mm_mul_ps(rinv20,rinv20);
 921             rinvsq30         = _mm_mul_ps(rinv30,rinv30);
 922
 923             /* Load parameters for j particles */
 924             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
 925                                                               charge+jnrC+0,charge+jnrD+0);
 926
 927             /**************************
 928              * CALCULATE INTERACTIONS *
 929              **************************/
 930
 931             r10              = _mm_mul_ps(rsq10,rinv10);
 932             r10              = _mm_andnot_ps(dummy_mask,r10);
 933
 934             /* Compute parameters for interactions between i and j atoms */
 935             qq10             = _mm_mul_ps(iq1,jq0);
 936
 937             /* EWALD ELECTROSTATICS */
 938
 939             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 940             ewrt             = _mm_mul_ps(r10,ewtabscale);
 941             ewitab           = _mm_cvttps_epi32(ewrt);
 942             eweps            = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
 943             gmx_mm_load_4pair_swizzle_ps(ewtab + gmx_mm_extract_epi32(ewitab,0),ewtab + gmx_mm_extract_epi32(ewitab,1),
 944                                          ewtab + gmx_mm_extract_epi32(ewitab,2),ewtab + gmx_mm_extract_epi32(ewitab,3),
 945                                          &ewtabF,&ewtabFn);
 946             felec            = _mm_add_ps(_mm_mul_ps( _mm_sub_ps(one,eweps),ewtabF),_mm_mul_ps(eweps,ewtabFn));
 947             felec            = _mm_mul_ps(_mm_mul_ps(qq10,rinv10),_mm_sub_ps(rinvsq10,felec));
 948
 949             fscal            = felec;
 950
 951             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 952
 953             /* Calculate temporary vectorial force */
 954             tx               = _mm_mul_ps(fscal,dx10);
 955             ty               = _mm_mul_ps(fscal,dy10);
 956             tz               = _mm_mul_ps(fscal,dz10);
 957
 958             /* Update vectorial force */
 959             fix1             = _mm_add_ps(fix1,tx);
 960             fiy1             = _mm_add_ps(fiy1,ty);
 961             fiz1             = _mm_add_ps(fiz1,tz);
 962
 963             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
 964             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
 965             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
 966             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
 967             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
 968
 969             /**************************
 970              * CALCULATE INTERACTIONS *
 971              **************************/
 972
 973             r20              = _mm_mul_ps(rsq20,rinv20);
 974             r20              = _mm_andnot_ps(dummy_mask,r20);
 975
 976             /* Compute parameters for interactions between i and j atoms */
 977             qq20             = _mm_mul_ps(iq2,jq0);
 978
 979             /* EWALD ELECTROSTATICS */
 980
 981             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 982             ewrt             = _mm_mul_ps(r20,ewtabscale);
 983             ewitab           = _mm_cvttps_epi32(ewrt);
 984             eweps            = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
 985             gmx_mm_load_4pair_swizzle_ps(ewtab + gmx_mm_extract_epi32(ewitab,0),ewtab + gmx_mm_extract_epi32(ewitab,1),
 986                                          ewtab + gmx_mm_extract_epi32(ewitab,2),ewtab + gmx_mm_extract_epi32(ewitab,3),
 987                                          &ewtabF,&ewtabFn);
 988             felec            = _mm_add_ps(_mm_mul_ps( _mm_sub_ps(one,eweps),ewtabF),_mm_mul_ps(eweps,ewtabFn));
 989             felec            = _mm_mul_ps(_mm_mul_ps(qq20,rinv20),_mm_sub_ps(rinvsq20,felec));
 990
 991             fscal            = felec;
 992
 993             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 994
 995             /* Calculate temporary vectorial force */
 996             tx               = _mm_mul_ps(fscal,dx20);
 997             ty               = _mm_mul_ps(fscal,dy20);
 998             tz               = _mm_mul_ps(fscal,dz20);
 999
1000             /* Update vectorial force */
1001             fix2             = _mm_add_ps(fix2,tx);
1002             fiy2             = _mm_add_ps(fiy2,ty);
1003             fiz2             = _mm_add_ps(fiz2,tz);
1004
1005             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
1006             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
1007             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
1008             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
1009             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
1010
1011             /**************************
1012              * CALCULATE INTERACTIONS *
1013              **************************/
1014
1015             r30              = _mm_mul_ps(rsq30,rinv30);
1016             r30              = _mm_andnot_ps(dummy_mask,r30);
1017
1018             /* Compute parameters for interactions between i and j atoms */
1019             qq30             = _mm_mul_ps(iq3,jq0);
1020
1021             /* EWALD ELECTROSTATICS */
1022
1023             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1024             ewrt             = _mm_mul_ps(r30,ewtabscale);
1025             ewitab           = _mm_cvttps_epi32(ewrt);
1026             eweps            = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
1027             gmx_mm_load_4pair_swizzle_ps(ewtab + gmx_mm_extract_epi32(ewitab,0),ewtab + gmx_mm_extract_epi32(ewitab,1),
1028                                          ewtab + gmx_mm_extract_epi32(ewitab,2),ewtab + gmx_mm_extract_epi32(ewitab,3),
1029                                          &ewtabF,&ewtabFn);
1030             felec            = _mm_add_ps(_mm_mul_ps( _mm_sub_ps(one,eweps),ewtabF),_mm_mul_ps(eweps,ewtabFn));
1031             felec            = _mm_mul_ps(_mm_mul_ps(qq30,rinv30),_mm_sub_ps(rinvsq30,felec));
1032
1033             fscal            = felec;
1034
1035             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1036
1037             /* Calculate temporary vectorial force */
1038             tx               = _mm_mul_ps(fscal,dx30);
1039             ty               = _mm_mul_ps(fscal,dy30);
1040             tz               = _mm_mul_ps(fscal,dz30);
1041
1042             /* Update vectorial force */
1043             fix3             = _mm_add_ps(fix3,tx);
1044             fiy3             = _mm_add_ps(fiy3,ty);
1045             fiz3             = _mm_add_ps(fiz3,tz);
1046
1047             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
1048             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
1049             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
1050             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
1051             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
1052
1053             /* Inner loop uses 111 flops */
1054         }
1055
1056         /* End of innermost loop */
1057
1058         gmx_mm_update_iforce_3atom_swizzle_ps(fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
1059                                               f+i_coord_offset+DIM,fshift+i_shift_offset);
1060
1061         /* Increment number of inner iterations */
1062         inneriter                  += j_index_end - j_index_start;
1063
1064         /* Outer loop uses 18 flops */
1065     }
1066
1067     /* Increment number of outer iterations */
1068     outeriter        += nri;
1069
1070     /* Update outer/inner flops */
1071
1072     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_W4_F,outeriter*18 + inneriter*111);
1073 }