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

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