src/gmxlib/nonbonded/nb_kernel_avx_256_double/nb_kernel_ElecEwSh_VdwNone_GeomW4W4_avx_256_double.c

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