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

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