src/gromacs/gmxlib/nonbonded/nb_kernel_avx_256_double/nb_kernel_ElecEwSh_VdwNone_GeomW4W4_avx_256_double.c

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