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

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