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

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*
  36  * Note: this file was generated by the GROMACS avx_256_double kernel generator.
  37  */
  38 #include "config.h"
  39
  40 #include <math.h>
  41
  42 #include "../nb_kernel.h"
  43 #include "types/simple.h"
  44 #include "gromacs/math/vec.h"
  45 #include "nrnb.h"
  46
  47 #include "gromacs/simd/math_x86_avx_256_double.h"
  48 #include "kernelutil_x86_avx_256_double.h"
  49
  50 /*
  51  * Gromacs nonbonded kernel:   nb_kernel_ElecCoul_VdwLJ_GeomW3W3_VF_avx_256_double
  52  * Electrostatics interaction: Coulomb
  53  * VdW interaction:            LennardJones
  54  * Geometry:                   Water3-Water3
  55  * Calculate force/pot:        PotentialAndForce
  56  */
  57 void
  58 nb_kernel_ElecCoul_VdwLJ_GeomW3W3_VF_avx_256_double
  59                     (t_nblist                    * gmx_restrict       nlist,
  60                      rvec                        * gmx_restrict          xx,
  61                      rvec                        * gmx_restrict          ff,
  62                      t_forcerec                  * gmx_restrict          fr,
  63                      t_mdatoms                   * gmx_restrict     mdatoms,
  64                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
  65                      t_nrnb                      * gmx_restrict        nrnb)
  66 {
  67     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
  68      * just 0 for non-waters.
  69      * Suffixes A,B,C,D refer to j loop unrolling done with AVX, e.g. for the four different
  70      * jnr indices corresponding to data put in the four positions in the SIMD register.
  71      */
  72     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
  73     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
  74     int              jnrA,jnrB,jnrC,jnrD;
  75     int              jnrlistA,jnrlistB,jnrlistC,jnrlistD;
  76     int              jnrlistE,jnrlistF,jnrlistG,jnrlistH;
  77     int              j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
  78     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
  79     real             rcutoff_scalar;
  80     real             *shiftvec,*fshift,*x,*f;
  81     real             *fjptrA,*fjptrB,*fjptrC,*fjptrD;
  82     real             scratch[4*DIM];
  83     __m256d          tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
  84     real *           vdwioffsetptr0;
  85     __m256d          ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
  86     real *           vdwioffsetptr1;
  87     __m256d          ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
  88     real *           vdwioffsetptr2;
  89     __m256d          ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
  90     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
  91     __m256d          jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
  92     int              vdwjidx1A,vdwjidx1B,vdwjidx1C,vdwjidx1D;
  93     __m256d          jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
  94     int              vdwjidx2A,vdwjidx2B,vdwjidx2C,vdwjidx2D;
  95     __m256d          jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
  96     __m256d          dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
  97     __m256d          dx01,dy01,dz01,rsq01,rinv01,rinvsq01,r01,qq01,c6_01,c12_01;
  98     __m256d          dx02,dy02,dz02,rsq02,rinv02,rinvsq02,r02,qq02,c6_02,c12_02;
  99     __m256d          dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
 100     __m256d          dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11;
 101     __m256d          dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12;
 102     __m256d          dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
 103     __m256d          dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21;
 104     __m256d          dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22;
 105     __m256d          velec,felec,velecsum,facel,crf,krf,krf2;
 106     real             *charge;
 107     int              nvdwtype;
 108     __m256d          rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
 109     int              *vdwtype;
 110     real             *vdwparam;
 111     __m256d          one_sixth   = _mm256_set1_pd(1.0/6.0);
 112     __m256d          one_twelfth = _mm256_set1_pd(1.0/12.0);
 113     __m256d          dummy_mask,cutoff_mask;
 114     __m128           tmpmask0,tmpmask1;
 115     __m256d          signbit = _mm256_castsi256_pd( _mm256_set1_epi32(0x80000000) );
 116     __m256d          one     = _mm256_set1_pd(1.0);
 117     __m256d          two     = _mm256_set1_pd(2.0);
 118     x                = xx[0];
 119     f                = ff[0];
 120
 121     nri              = nlist->nri;
 122     iinr             = nlist->iinr;
 123     jindex           = nlist->jindex;
 124     jjnr             = nlist->jjnr;
 125     shiftidx         = nlist->shift;
 126     gid              = nlist->gid;
 127     shiftvec         = fr->shift_vec[0];
 128     fshift           = fr->fshift[0];
 129     facel            = _mm256_set1_pd(fr->epsfac);
 130     charge           = mdatoms->chargeA;
 131     nvdwtype         = fr->ntype;
 132     vdwparam         = fr->nbfp;
 133     vdwtype          = mdatoms->typeA;
 134
 135     /* Setup water-specific parameters */
 136     inr              = nlist->iinr[0];
 137     iq0              = _mm256_mul_pd(facel,_mm256_set1_pd(charge[inr+0]));
 138     iq1              = _mm256_mul_pd(facel,_mm256_set1_pd(charge[inr+1]));
 139     iq2              = _mm256_mul_pd(facel,_mm256_set1_pd(charge[inr+2]));
 140     vdwioffsetptr0   = vdwparam+2*nvdwtype*vdwtype[inr+0];
 141
 142     jq0              = _mm256_set1_pd(charge[inr+0]);
 143     jq1              = _mm256_set1_pd(charge[inr+1]);
 144     jq2              = _mm256_set1_pd(charge[inr+2]);
 145     vdwjidx0A        = 2*vdwtype[inr+0];
 146     qq00             = _mm256_mul_pd(iq0,jq0);
 147     c6_00            = _mm256_set1_pd(vdwioffsetptr0[vdwjidx0A]);
 148     c12_00           = _mm256_set1_pd(vdwioffsetptr0[vdwjidx0A+1]);
 149     qq01             = _mm256_mul_pd(iq0,jq1);
 150     qq02             = _mm256_mul_pd(iq0,jq2);
 151     qq10             = _mm256_mul_pd(iq1,jq0);
 152     qq11             = _mm256_mul_pd(iq1,jq1);
 153     qq12             = _mm256_mul_pd(iq1,jq2);
 154     qq20             = _mm256_mul_pd(iq2,jq0);
 155     qq21             = _mm256_mul_pd(iq2,jq1);
 156     qq22             = _mm256_mul_pd(iq2,jq2);
 157
 158     /* Avoid stupid compiler warnings */
 159     jnrA = jnrB = jnrC = jnrD = 0;
 160     j_coord_offsetA = 0;
 161     j_coord_offsetB = 0;
 162     j_coord_offsetC = 0;
 163     j_coord_offsetD = 0;
 164
 165     outeriter        = 0;
 166     inneriter        = 0;
 167
 168     for(iidx=0;iidx<4*DIM;iidx++)
 169     {
 170         scratch[iidx] = 0.0;
 171     }
 172
 173     /* Start outer loop over neighborlists */
 174     for(iidx=0; iidx<nri; iidx++)
 175     {
 176         /* Load shift vector for this list */
 177         i_shift_offset   = DIM*shiftidx[iidx];
 178
 179         /* Load limits for loop over neighbors */
 180         j_index_start    = jindex[iidx];
 181         j_index_end      = jindex[iidx+1];
 182
 183         /* Get outer coordinate index */
 184         inr              = iinr[iidx];
 185         i_coord_offset   = DIM*inr;
 186
 187         /* Load i particle coords and add shift vector */
 188         gmx_mm256_load_shift_and_3rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,
 189                                                     &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2);
 190
 191         fix0             = _mm256_setzero_pd();
 192         fiy0             = _mm256_setzero_pd();
 193         fiz0             = _mm256_setzero_pd();
 194         fix1             = _mm256_setzero_pd();
 195         fiy1             = _mm256_setzero_pd();
 196         fiz1             = _mm256_setzero_pd();
 197         fix2             = _mm256_setzero_pd();
 198         fiy2             = _mm256_setzero_pd();
 199         fiz2             = _mm256_setzero_pd();
 200
 201         /* Reset potential sums */
 202         velecsum         = _mm256_setzero_pd();
 203         vvdwsum          = _mm256_setzero_pd();
 204
 205         /* Start inner kernel loop */
 206         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
 207         {
 208
 209             /* Get j neighbor index, and coordinate index */
 210             jnrA             = jjnr[jidx];
 211             jnrB             = jjnr[jidx+1];
 212             jnrC             = jjnr[jidx+2];
 213             jnrD             = jjnr[jidx+3];
 214             j_coord_offsetA  = DIM*jnrA;
 215             j_coord_offsetB  = DIM*jnrB;
 216             j_coord_offsetC  = DIM*jnrC;
 217             j_coord_offsetD  = DIM*jnrD;
 218
 219             /* load j atom coordinates */
 220             gmx_mm256_load_3rvec_4ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
 221                                                  x+j_coord_offsetC,x+j_coord_offsetD,
 222                                               &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,&jy2,&jz2);
 223
 224             /* Calculate displacement vector */
 225             dx00             = _mm256_sub_pd(ix0,jx0);
 226             dy00             = _mm256_sub_pd(iy0,jy0);
 227             dz00             = _mm256_sub_pd(iz0,jz0);
 228             dx01             = _mm256_sub_pd(ix0,jx1);
 229             dy01             = _mm256_sub_pd(iy0,jy1);
 230             dz01             = _mm256_sub_pd(iz0,jz1);
 231             dx02             = _mm256_sub_pd(ix0,jx2);
 232             dy02             = _mm256_sub_pd(iy0,jy2);
 233             dz02             = _mm256_sub_pd(iz0,jz2);
 234             dx10             = _mm256_sub_pd(ix1,jx0);
 235             dy10             = _mm256_sub_pd(iy1,jy0);
 236             dz10             = _mm256_sub_pd(iz1,jz0);
 237             dx11             = _mm256_sub_pd(ix1,jx1);
 238             dy11             = _mm256_sub_pd(iy1,jy1);
 239             dz11             = _mm256_sub_pd(iz1,jz1);
 240             dx12             = _mm256_sub_pd(ix1,jx2);
 241             dy12             = _mm256_sub_pd(iy1,jy2);
 242             dz12             = _mm256_sub_pd(iz1,jz2);
 243             dx20             = _mm256_sub_pd(ix2,jx0);
 244             dy20             = _mm256_sub_pd(iy2,jy0);
 245             dz20             = _mm256_sub_pd(iz2,jz0);
 246             dx21             = _mm256_sub_pd(ix2,jx1);
 247             dy21             = _mm256_sub_pd(iy2,jy1);
 248             dz21             = _mm256_sub_pd(iz2,jz1);
 249             dx22             = _mm256_sub_pd(ix2,jx2);
 250             dy22             = _mm256_sub_pd(iy2,jy2);
 251             dz22             = _mm256_sub_pd(iz2,jz2);
 252
 253             /* Calculate squared distance and things based on it */
 254             rsq00            = gmx_mm256_calc_rsq_pd(dx00,dy00,dz00);
 255             rsq01            = gmx_mm256_calc_rsq_pd(dx01,dy01,dz01);
 256             rsq02            = gmx_mm256_calc_rsq_pd(dx02,dy02,dz02);
 257             rsq10            = gmx_mm256_calc_rsq_pd(dx10,dy10,dz10);
 258             rsq11            = gmx_mm256_calc_rsq_pd(dx11,dy11,dz11);
 259             rsq12            = gmx_mm256_calc_rsq_pd(dx12,dy12,dz12);
 260             rsq20            = gmx_mm256_calc_rsq_pd(dx20,dy20,dz20);
 261             rsq21            = gmx_mm256_calc_rsq_pd(dx21,dy21,dz21);
 262             rsq22            = gmx_mm256_calc_rsq_pd(dx22,dy22,dz22);
 263
 264             rinv00           = gmx_mm256_invsqrt_pd(rsq00);
 265             rinv01           = gmx_mm256_invsqrt_pd(rsq01);
 266             rinv02           = gmx_mm256_invsqrt_pd(rsq02);
 267             rinv10           = gmx_mm256_invsqrt_pd(rsq10);
 268             rinv11           = gmx_mm256_invsqrt_pd(rsq11);
 269             rinv12           = gmx_mm256_invsqrt_pd(rsq12);
 270             rinv20           = gmx_mm256_invsqrt_pd(rsq20);
 271             rinv21           = gmx_mm256_invsqrt_pd(rsq21);
 272             rinv22           = gmx_mm256_invsqrt_pd(rsq22);
 273
 274             rinvsq00         = _mm256_mul_pd(rinv00,rinv00);
 275             rinvsq01         = _mm256_mul_pd(rinv01,rinv01);
 276             rinvsq02         = _mm256_mul_pd(rinv02,rinv02);
 277             rinvsq10         = _mm256_mul_pd(rinv10,rinv10);
 278             rinvsq11         = _mm256_mul_pd(rinv11,rinv11);
 279             rinvsq12         = _mm256_mul_pd(rinv12,rinv12);
 280             rinvsq20         = _mm256_mul_pd(rinv20,rinv20);
 281             rinvsq21         = _mm256_mul_pd(rinv21,rinv21);
 282             rinvsq22         = _mm256_mul_pd(rinv22,rinv22);
 283
 284             fjx0             = _mm256_setzero_pd();
 285             fjy0             = _mm256_setzero_pd();
 286             fjz0             = _mm256_setzero_pd();
 287             fjx1             = _mm256_setzero_pd();
 288             fjy1             = _mm256_setzero_pd();
 289             fjz1             = _mm256_setzero_pd();
 290             fjx2             = _mm256_setzero_pd();
 291             fjy2             = _mm256_setzero_pd();
 292             fjz2             = _mm256_setzero_pd();
 293
 294             /**************************
 295              * CALCULATE INTERACTIONS *
 296              **************************/
 297
 298             /* COULOMB ELECTROSTATICS */
 299             velec            = _mm256_mul_pd(qq00,rinv00);
 300             felec            = _mm256_mul_pd(velec,rinvsq00);
 301
 302             /* LENNARD-JONES DISPERSION/REPULSION */
 303
 304             rinvsix          = _mm256_mul_pd(_mm256_mul_pd(rinvsq00,rinvsq00),rinvsq00);
 305             vvdw6            = _mm256_mul_pd(c6_00,rinvsix);
 306             vvdw12           = _mm256_mul_pd(c12_00,_mm256_mul_pd(rinvsix,rinvsix));
 307             vvdw             = _mm256_sub_pd( _mm256_mul_pd(vvdw12,one_twelfth) , _mm256_mul_pd(vvdw6,one_sixth) );
 308             fvdw             = _mm256_mul_pd(_mm256_sub_pd(vvdw12,vvdw6),rinvsq00);
 309
 310             /* Update potential sum for this i atom from the interaction with this j atom. */
 311             velecsum         = _mm256_add_pd(velecsum,velec);
 312             vvdwsum          = _mm256_add_pd(vvdwsum,vvdw);
 313
 314             fscal            = _mm256_add_pd(felec,fvdw);
 315
 316             /* Calculate temporary vectorial force */
 317             tx               = _mm256_mul_pd(fscal,dx00);
 318             ty               = _mm256_mul_pd(fscal,dy00);
 319             tz               = _mm256_mul_pd(fscal,dz00);
 320
 321             /* Update vectorial force */
 322             fix0             = _mm256_add_pd(fix0,tx);
 323             fiy0             = _mm256_add_pd(fiy0,ty);
 324             fiz0             = _mm256_add_pd(fiz0,tz);
 325
 326             fjx0             = _mm256_add_pd(fjx0,tx);
 327             fjy0             = _mm256_add_pd(fjy0,ty);
 328             fjz0             = _mm256_add_pd(fjz0,tz);
 329
 330             /**************************
 331              * CALCULATE INTERACTIONS *
 332              **************************/
 333
 334             /* COULOMB ELECTROSTATICS */
 335             velec            = _mm256_mul_pd(qq01,rinv01);
 336             felec            = _mm256_mul_pd(velec,rinvsq01);
 337
 338             /* Update potential sum for this i atom from the interaction with this j atom. */
 339             velecsum         = _mm256_add_pd(velecsum,velec);
 340
 341             fscal            = felec;
 342
 343             /* Calculate temporary vectorial force */
 344             tx               = _mm256_mul_pd(fscal,dx01);
 345             ty               = _mm256_mul_pd(fscal,dy01);
 346             tz               = _mm256_mul_pd(fscal,dz01);
 347
 348             /* Update vectorial force */
 349             fix0             = _mm256_add_pd(fix0,tx);
 350             fiy0             = _mm256_add_pd(fiy0,ty);
 351             fiz0             = _mm256_add_pd(fiz0,tz);
 352
 353             fjx1             = _mm256_add_pd(fjx1,tx);
 354             fjy1             = _mm256_add_pd(fjy1,ty);
 355             fjz1             = _mm256_add_pd(fjz1,tz);
 356
 357             /**************************
 358              * CALCULATE INTERACTIONS *
 359              **************************/
 360
 361             /* COULOMB ELECTROSTATICS */
 362             velec            = _mm256_mul_pd(qq02,rinv02);
 363             felec            = _mm256_mul_pd(velec,rinvsq02);
 364
 365             /* Update potential sum for this i atom from the interaction with this j atom. */
 366             velecsum         = _mm256_add_pd(velecsum,velec);
 367
 368             fscal            = felec;
 369
 370             /* Calculate temporary vectorial force */
 371             tx               = _mm256_mul_pd(fscal,dx02);
 372             ty               = _mm256_mul_pd(fscal,dy02);
 373             tz               = _mm256_mul_pd(fscal,dz02);
 374
 375             /* Update vectorial force */
 376             fix0             = _mm256_add_pd(fix0,tx);
 377             fiy0             = _mm256_add_pd(fiy0,ty);
 378             fiz0             = _mm256_add_pd(fiz0,tz);
 379
 380             fjx2             = _mm256_add_pd(fjx2,tx);
 381             fjy2             = _mm256_add_pd(fjy2,ty);
 382             fjz2             = _mm256_add_pd(fjz2,tz);
 383
 384             /**************************
 385              * CALCULATE INTERACTIONS *
 386              **************************/
 387
 388             /* COULOMB ELECTROSTATICS */
 389             velec            = _mm256_mul_pd(qq10,rinv10);
 390             felec            = _mm256_mul_pd(velec,rinvsq10);
 391
 392             /* Update potential sum for this i atom from the interaction with this j atom. */
 393             velecsum         = _mm256_add_pd(velecsum,velec);
 394
 395             fscal            = felec;
 396
 397             /* Calculate temporary vectorial force */
 398             tx               = _mm256_mul_pd(fscal,dx10);
 399             ty               = _mm256_mul_pd(fscal,dy10);
 400             tz               = _mm256_mul_pd(fscal,dz10);
 401
 402             /* Update vectorial force */
 403             fix1             = _mm256_add_pd(fix1,tx);
 404             fiy1             = _mm256_add_pd(fiy1,ty);
 405             fiz1             = _mm256_add_pd(fiz1,tz);
 406
 407             fjx0             = _mm256_add_pd(fjx0,tx);
 408             fjy0             = _mm256_add_pd(fjy0,ty);
 409             fjz0             = _mm256_add_pd(fjz0,tz);
 410
 411             /**************************
 412              * CALCULATE INTERACTIONS *
 413              **************************/
 414
 415             /* COULOMB ELECTROSTATICS */
 416             velec            = _mm256_mul_pd(qq11,rinv11);
 417             felec            = _mm256_mul_pd(velec,rinvsq11);
 418
 419             /* Update potential sum for this i atom from the interaction with this j atom. */
 420             velecsum         = _mm256_add_pd(velecsum,velec);
 421
 422             fscal            = felec;
 423
 424             /* Calculate temporary vectorial force */
 425             tx               = _mm256_mul_pd(fscal,dx11);
 426             ty               = _mm256_mul_pd(fscal,dy11);
 427             tz               = _mm256_mul_pd(fscal,dz11);
 428
 429             /* Update vectorial force */
 430             fix1             = _mm256_add_pd(fix1,tx);
 431             fiy1             = _mm256_add_pd(fiy1,ty);
 432             fiz1             = _mm256_add_pd(fiz1,tz);
 433
 434             fjx1             = _mm256_add_pd(fjx1,tx);
 435             fjy1             = _mm256_add_pd(fjy1,ty);
 436             fjz1             = _mm256_add_pd(fjz1,tz);
 437
 438             /**************************
 439              * CALCULATE INTERACTIONS *
 440              **************************/
 441
 442             /* COULOMB ELECTROSTATICS */
 443             velec            = _mm256_mul_pd(qq12,rinv12);
 444             felec            = _mm256_mul_pd(velec,rinvsq12);
 445
 446             /* Update potential sum for this i atom from the interaction with this j atom. */
 447             velecsum         = _mm256_add_pd(velecsum,velec);
 448
 449             fscal            = felec;
 450
 451             /* Calculate temporary vectorial force */
 452             tx               = _mm256_mul_pd(fscal,dx12);
 453             ty               = _mm256_mul_pd(fscal,dy12);
 454             tz               = _mm256_mul_pd(fscal,dz12);
 455
 456             /* Update vectorial force */
 457             fix1             = _mm256_add_pd(fix1,tx);
 458             fiy1             = _mm256_add_pd(fiy1,ty);
 459             fiz1             = _mm256_add_pd(fiz1,tz);
 460
 461             fjx2             = _mm256_add_pd(fjx2,tx);
 462             fjy2             = _mm256_add_pd(fjy2,ty);
 463             fjz2             = _mm256_add_pd(fjz2,tz);
 464
 465             /**************************
 466              * CALCULATE INTERACTIONS *
 467              **************************/
 468
 469             /* COULOMB ELECTROSTATICS */
 470             velec            = _mm256_mul_pd(qq20,rinv20);
 471             felec            = _mm256_mul_pd(velec,rinvsq20);
 472
 473             /* Update potential sum for this i atom from the interaction with this j atom. */
 474             velecsum         = _mm256_add_pd(velecsum,velec);
 475
 476             fscal            = felec;
 477
 478             /* Calculate temporary vectorial force */
 479             tx               = _mm256_mul_pd(fscal,dx20);
 480             ty               = _mm256_mul_pd(fscal,dy20);
 481             tz               = _mm256_mul_pd(fscal,dz20);
 482
 483             /* Update vectorial force */
 484             fix2             = _mm256_add_pd(fix2,tx);
 485             fiy2             = _mm256_add_pd(fiy2,ty);
 486             fiz2             = _mm256_add_pd(fiz2,tz);
 487
 488             fjx0             = _mm256_add_pd(fjx0,tx);
 489             fjy0             = _mm256_add_pd(fjy0,ty);
 490             fjz0             = _mm256_add_pd(fjz0,tz);
 491
 492             /**************************
 493              * CALCULATE INTERACTIONS *
 494              **************************/
 495
 496             /* COULOMB ELECTROSTATICS */
 497             velec            = _mm256_mul_pd(qq21,rinv21);
 498             felec            = _mm256_mul_pd(velec,rinvsq21);
 499
 500             /* Update potential sum for this i atom from the interaction with this j atom. */
 501             velecsum         = _mm256_add_pd(velecsum,velec);
 502
 503             fscal            = felec;
 504
 505             /* Calculate temporary vectorial force */
 506             tx               = _mm256_mul_pd(fscal,dx21);
 507             ty               = _mm256_mul_pd(fscal,dy21);
 508             tz               = _mm256_mul_pd(fscal,dz21);
 509
 510             /* Update vectorial force */
 511             fix2             = _mm256_add_pd(fix2,tx);
 512             fiy2             = _mm256_add_pd(fiy2,ty);
 513             fiz2             = _mm256_add_pd(fiz2,tz);
 514
 515             fjx1             = _mm256_add_pd(fjx1,tx);
 516             fjy1             = _mm256_add_pd(fjy1,ty);
 517             fjz1             = _mm256_add_pd(fjz1,tz);
 518
 519             /**************************
 520              * CALCULATE INTERACTIONS *
 521              **************************/
 522
 523             /* COULOMB ELECTROSTATICS */
 524             velec            = _mm256_mul_pd(qq22,rinv22);
 525             felec            = _mm256_mul_pd(velec,rinvsq22);
 526
 527             /* Update potential sum for this i atom from the interaction with this j atom. */
 528             velecsum         = _mm256_add_pd(velecsum,velec);
 529
 530             fscal            = felec;
 531
 532             /* Calculate temporary vectorial force */
 533             tx               = _mm256_mul_pd(fscal,dx22);
 534             ty               = _mm256_mul_pd(fscal,dy22);
 535             tz               = _mm256_mul_pd(fscal,dz22);
 536
 537             /* Update vectorial force */
 538             fix2             = _mm256_add_pd(fix2,tx);
 539             fiy2             = _mm256_add_pd(fiy2,ty);
 540             fiz2             = _mm256_add_pd(fiz2,tz);
 541
 542             fjx2             = _mm256_add_pd(fjx2,tx);
 543             fjy2             = _mm256_add_pd(fjy2,ty);
 544             fjz2             = _mm256_add_pd(fjz2,tz);
 545
 546             fjptrA             = f+j_coord_offsetA;
 547             fjptrB             = f+j_coord_offsetB;
 548             fjptrC             = f+j_coord_offsetC;
 549             fjptrD             = f+j_coord_offsetD;
 550
 551             gmx_mm256_decrement_3rvec_4ptr_swizzle_pd(fjptrA,fjptrB,fjptrC,fjptrD,
 552                                                       fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,fjx2,fjy2,fjz2);
 553
 554             /* Inner loop uses 255 flops */
 555         }
 556
 557         if(jidx<j_index_end)
 558         {
 559
 560             /* Get j neighbor index, and coordinate index */
 561             jnrlistA         = jjnr[jidx];
 562             jnrlistB         = jjnr[jidx+1];
 563             jnrlistC         = jjnr[jidx+2];
 564             jnrlistD         = jjnr[jidx+3];
 565             /* Sign of each element will be negative for non-real atoms.
 566              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
 567              * so use it as val = _mm_andnot_pd(mask,val) to clear dummy entries.
 568              */
 569             tmpmask0 = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
 570
 571             tmpmask1 = _mm_permute_ps(tmpmask0,_GMX_MM_PERMUTE(3,3,2,2));
 572             tmpmask0 = _mm_permute_ps(tmpmask0,_GMX_MM_PERMUTE(1,1,0,0));
 573             dummy_mask = _mm256_castps_pd(gmx_mm256_set_m128(tmpmask1,tmpmask0));
 574
 575             jnrA       = (jnrlistA>=0) ? jnrlistA : 0;
 576             jnrB       = (jnrlistB>=0) ? jnrlistB : 0;
 577             jnrC       = (jnrlistC>=0) ? jnrlistC : 0;
 578             jnrD       = (jnrlistD>=0) ? jnrlistD : 0;
 579             j_coord_offsetA  = DIM*jnrA;
 580             j_coord_offsetB  = DIM*jnrB;
 581             j_coord_offsetC  = DIM*jnrC;
 582             j_coord_offsetD  = DIM*jnrD;
 583
 584             /* load j atom coordinates */
 585             gmx_mm256_load_3rvec_4ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
 586                                                  x+j_coord_offsetC,x+j_coord_offsetD,
 587                                               &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,&jy2,&jz2);
 588
 589             /* Calculate displacement vector */
 590             dx00             = _mm256_sub_pd(ix0,jx0);
 591             dy00             = _mm256_sub_pd(iy0,jy0);
 592             dz00             = _mm256_sub_pd(iz0,jz0);
 593             dx01             = _mm256_sub_pd(ix0,jx1);
 594             dy01             = _mm256_sub_pd(iy0,jy1);
 595             dz01             = _mm256_sub_pd(iz0,jz1);
 596             dx02             = _mm256_sub_pd(ix0,jx2);
 597             dy02             = _mm256_sub_pd(iy0,jy2);
 598             dz02             = _mm256_sub_pd(iz0,jz2);
 599             dx10             = _mm256_sub_pd(ix1,jx0);
 600             dy10             = _mm256_sub_pd(iy1,jy0);
 601             dz10             = _mm256_sub_pd(iz1,jz0);
 602             dx11             = _mm256_sub_pd(ix1,jx1);
 603             dy11             = _mm256_sub_pd(iy1,jy1);
 604             dz11             = _mm256_sub_pd(iz1,jz1);
 605             dx12             = _mm256_sub_pd(ix1,jx2);
 606             dy12             = _mm256_sub_pd(iy1,jy2);
 607             dz12             = _mm256_sub_pd(iz1,jz2);
 608             dx20             = _mm256_sub_pd(ix2,jx0);
 609             dy20             = _mm256_sub_pd(iy2,jy0);
 610             dz20             = _mm256_sub_pd(iz2,jz0);
 611             dx21             = _mm256_sub_pd(ix2,jx1);
 612             dy21             = _mm256_sub_pd(iy2,jy1);
 613             dz21             = _mm256_sub_pd(iz2,jz1);
 614             dx22             = _mm256_sub_pd(ix2,jx2);
 615             dy22             = _mm256_sub_pd(iy2,jy2);
 616             dz22             = _mm256_sub_pd(iz2,jz2);
 617
 618             /* Calculate squared distance and things based on it */
 619             rsq00            = gmx_mm256_calc_rsq_pd(dx00,dy00,dz00);
 620             rsq01            = gmx_mm256_calc_rsq_pd(dx01,dy01,dz01);
 621             rsq02            = gmx_mm256_calc_rsq_pd(dx02,dy02,dz02);
 622             rsq10            = gmx_mm256_calc_rsq_pd(dx10,dy10,dz10);
 623             rsq11            = gmx_mm256_calc_rsq_pd(dx11,dy11,dz11);
 624             rsq12            = gmx_mm256_calc_rsq_pd(dx12,dy12,dz12);
 625             rsq20            = gmx_mm256_calc_rsq_pd(dx20,dy20,dz20);
 626             rsq21            = gmx_mm256_calc_rsq_pd(dx21,dy21,dz21);
 627             rsq22            = gmx_mm256_calc_rsq_pd(dx22,dy22,dz22);
 628
 629             rinv00           = gmx_mm256_invsqrt_pd(rsq00);
 630             rinv01           = gmx_mm256_invsqrt_pd(rsq01);
 631             rinv02           = gmx_mm256_invsqrt_pd(rsq02);
 632             rinv10           = gmx_mm256_invsqrt_pd(rsq10);
 633             rinv11           = gmx_mm256_invsqrt_pd(rsq11);
 634             rinv12           = gmx_mm256_invsqrt_pd(rsq12);
 635             rinv20           = gmx_mm256_invsqrt_pd(rsq20);
 636             rinv21           = gmx_mm256_invsqrt_pd(rsq21);
 637             rinv22           = gmx_mm256_invsqrt_pd(rsq22);
 638
 639             rinvsq00         = _mm256_mul_pd(rinv00,rinv00);
 640             rinvsq01         = _mm256_mul_pd(rinv01,rinv01);
 641             rinvsq02         = _mm256_mul_pd(rinv02,rinv02);
 642             rinvsq10         = _mm256_mul_pd(rinv10,rinv10);
 643             rinvsq11         = _mm256_mul_pd(rinv11,rinv11);
 644             rinvsq12         = _mm256_mul_pd(rinv12,rinv12);
 645             rinvsq20         = _mm256_mul_pd(rinv20,rinv20);
 646             rinvsq21         = _mm256_mul_pd(rinv21,rinv21);
 647             rinvsq22         = _mm256_mul_pd(rinv22,rinv22);
 648
 649             fjx0             = _mm256_setzero_pd();
 650             fjy0             = _mm256_setzero_pd();
 651             fjz0             = _mm256_setzero_pd();
 652             fjx1             = _mm256_setzero_pd();
 653             fjy1             = _mm256_setzero_pd();
 654             fjz1             = _mm256_setzero_pd();
 655             fjx2             = _mm256_setzero_pd();
 656             fjy2             = _mm256_setzero_pd();
 657             fjz2             = _mm256_setzero_pd();
 658
 659             /**************************
 660              * CALCULATE INTERACTIONS *
 661              **************************/
 662
 663             /* COULOMB ELECTROSTATICS */
 664             velec            = _mm256_mul_pd(qq00,rinv00);
 665             felec            = _mm256_mul_pd(velec,rinvsq00);
 666
 667             /* LENNARD-JONES DISPERSION/REPULSION */
 668
 669             rinvsix          = _mm256_mul_pd(_mm256_mul_pd(rinvsq00,rinvsq00),rinvsq00);
 670             vvdw6            = _mm256_mul_pd(c6_00,rinvsix);
 671             vvdw12           = _mm256_mul_pd(c12_00,_mm256_mul_pd(rinvsix,rinvsix));
 672             vvdw             = _mm256_sub_pd( _mm256_mul_pd(vvdw12,one_twelfth) , _mm256_mul_pd(vvdw6,one_sixth) );
 673             fvdw             = _mm256_mul_pd(_mm256_sub_pd(vvdw12,vvdw6),rinvsq00);
 674
 675             /* Update potential sum for this i atom from the interaction with this j atom. */
 676             velec            = _mm256_andnot_pd(dummy_mask,velec);
 677             velecsum         = _mm256_add_pd(velecsum,velec);
 678             vvdw             = _mm256_andnot_pd(dummy_mask,vvdw);
 679             vvdwsum          = _mm256_add_pd(vvdwsum,vvdw);
 680
 681             fscal            = _mm256_add_pd(felec,fvdw);
 682
 683             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
 684
 685             /* Calculate temporary vectorial force */
 686             tx               = _mm256_mul_pd(fscal,dx00);
 687             ty               = _mm256_mul_pd(fscal,dy00);
 688             tz               = _mm256_mul_pd(fscal,dz00);
 689
 690             /* Update vectorial force */
 691             fix0             = _mm256_add_pd(fix0,tx);
 692             fiy0             = _mm256_add_pd(fiy0,ty);
 693             fiz0             = _mm256_add_pd(fiz0,tz);
 694
 695             fjx0             = _mm256_add_pd(fjx0,tx);
 696             fjy0             = _mm256_add_pd(fjy0,ty);
 697             fjz0             = _mm256_add_pd(fjz0,tz);
 698
 699             /**************************
 700              * CALCULATE INTERACTIONS *
 701              **************************/
 702
 703             /* COULOMB ELECTROSTATICS */
 704             velec            = _mm256_mul_pd(qq01,rinv01);
 705             felec            = _mm256_mul_pd(velec,rinvsq01);
 706
 707             /* Update potential sum for this i atom from the interaction with this j atom. */
 708             velec            = _mm256_andnot_pd(dummy_mask,velec);
 709             velecsum         = _mm256_add_pd(velecsum,velec);
 710
 711             fscal            = felec;
 712
 713             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
 714
 715             /* Calculate temporary vectorial force */
 716             tx               = _mm256_mul_pd(fscal,dx01);
 717             ty               = _mm256_mul_pd(fscal,dy01);
 718             tz               = _mm256_mul_pd(fscal,dz01);
 719
 720             /* Update vectorial force */
 721             fix0             = _mm256_add_pd(fix0,tx);
 722             fiy0             = _mm256_add_pd(fiy0,ty);
 723             fiz0             = _mm256_add_pd(fiz0,tz);
 724
 725             fjx1             = _mm256_add_pd(fjx1,tx);
 726             fjy1             = _mm256_add_pd(fjy1,ty);
 727             fjz1             = _mm256_add_pd(fjz1,tz);
 728
 729             /**************************
 730              * CALCULATE INTERACTIONS *
 731              **************************/
 732
 733             /* COULOMB ELECTROSTATICS */
 734             velec            = _mm256_mul_pd(qq02,rinv02);
 735             felec            = _mm256_mul_pd(velec,rinvsq02);
 736
 737             /* Update potential sum for this i atom from the interaction with this j atom. */
 738             velec            = _mm256_andnot_pd(dummy_mask,velec);
 739             velecsum         = _mm256_add_pd(velecsum,velec);
 740
 741             fscal            = felec;
 742
 743             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
 744
 745             /* Calculate temporary vectorial force */
 746             tx               = _mm256_mul_pd(fscal,dx02);
 747             ty               = _mm256_mul_pd(fscal,dy02);
 748             tz               = _mm256_mul_pd(fscal,dz02);
 749
 750             /* Update vectorial force */
 751             fix0             = _mm256_add_pd(fix0,tx);
 752             fiy0             = _mm256_add_pd(fiy0,ty);
 753             fiz0             = _mm256_add_pd(fiz0,tz);
 754
 755             fjx2             = _mm256_add_pd(fjx2,tx);
 756             fjy2             = _mm256_add_pd(fjy2,ty);
 757             fjz2             = _mm256_add_pd(fjz2,tz);
 758
 759             /**************************
 760              * CALCULATE INTERACTIONS *
 761              **************************/
 762
 763             /* COULOMB ELECTROSTATICS */
 764             velec            = _mm256_mul_pd(qq10,rinv10);
 765             felec            = _mm256_mul_pd(velec,rinvsq10);
 766
 767             /* Update potential sum for this i atom from the interaction with this j atom. */
 768             velec            = _mm256_andnot_pd(dummy_mask,velec);
 769             velecsum         = _mm256_add_pd(velecsum,velec);
 770
 771             fscal            = felec;
 772
 773             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
 774
 775             /* Calculate temporary vectorial force */
 776             tx               = _mm256_mul_pd(fscal,dx10);
 777             ty               = _mm256_mul_pd(fscal,dy10);
 778             tz               = _mm256_mul_pd(fscal,dz10);
 779
 780             /* Update vectorial force */
 781             fix1             = _mm256_add_pd(fix1,tx);
 782             fiy1             = _mm256_add_pd(fiy1,ty);
 783             fiz1             = _mm256_add_pd(fiz1,tz);
 784
 785             fjx0             = _mm256_add_pd(fjx0,tx);
 786             fjy0             = _mm256_add_pd(fjy0,ty);
 787             fjz0             = _mm256_add_pd(fjz0,tz);
 788
 789             /**************************
 790              * CALCULATE INTERACTIONS *
 791              **************************/
 792
 793             /* COULOMB ELECTROSTATICS */
 794             velec            = _mm256_mul_pd(qq11,rinv11);
 795             felec            = _mm256_mul_pd(velec,rinvsq11);
 796
 797             /* Update potential sum for this i atom from the interaction with this j atom. */
 798             velec            = _mm256_andnot_pd(dummy_mask,velec);
 799             velecsum         = _mm256_add_pd(velecsum,velec);
 800
 801             fscal            = felec;
 802
 803             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
 804
 805             /* Calculate temporary vectorial force */
 806             tx               = _mm256_mul_pd(fscal,dx11);
 807             ty               = _mm256_mul_pd(fscal,dy11);
 808             tz               = _mm256_mul_pd(fscal,dz11);
 809
 810             /* Update vectorial force */
 811             fix1             = _mm256_add_pd(fix1,tx);
 812             fiy1             = _mm256_add_pd(fiy1,ty);
 813             fiz1             = _mm256_add_pd(fiz1,tz);
 814
 815             fjx1             = _mm256_add_pd(fjx1,tx);
 816             fjy1             = _mm256_add_pd(fjy1,ty);
 817             fjz1             = _mm256_add_pd(fjz1,tz);
 818
 819             /**************************
 820              * CALCULATE INTERACTIONS *
 821              **************************/
 822
 823             /* COULOMB ELECTROSTATICS */
 824             velec            = _mm256_mul_pd(qq12,rinv12);
 825             felec            = _mm256_mul_pd(velec,rinvsq12);
 826
 827             /* Update potential sum for this i atom from the interaction with this j atom. */
 828             velec            = _mm256_andnot_pd(dummy_mask,velec);
 829             velecsum         = _mm256_add_pd(velecsum,velec);
 830
 831             fscal            = felec;
 832
 833             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
 834
 835             /* Calculate temporary vectorial force */
 836             tx               = _mm256_mul_pd(fscal,dx12);
 837             ty               = _mm256_mul_pd(fscal,dy12);
 838             tz               = _mm256_mul_pd(fscal,dz12);
 839
 840             /* Update vectorial force */
 841             fix1             = _mm256_add_pd(fix1,tx);
 842             fiy1             = _mm256_add_pd(fiy1,ty);
 843             fiz1             = _mm256_add_pd(fiz1,tz);
 844
 845             fjx2             = _mm256_add_pd(fjx2,tx);
 846             fjy2             = _mm256_add_pd(fjy2,ty);
 847             fjz2             = _mm256_add_pd(fjz2,tz);
 848
 849             /**************************
 850              * CALCULATE INTERACTIONS *
 851              **************************/
 852
 853             /* COULOMB ELECTROSTATICS */
 854             velec            = _mm256_mul_pd(qq20,rinv20);
 855             felec            = _mm256_mul_pd(velec,rinvsq20);
 856
 857             /* Update potential sum for this i atom from the interaction with this j atom. */
 858             velec            = _mm256_andnot_pd(dummy_mask,velec);
 859             velecsum         = _mm256_add_pd(velecsum,velec);
 860
 861             fscal            = felec;
 862
 863             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
 864
 865             /* Calculate temporary vectorial force */
 866             tx               = _mm256_mul_pd(fscal,dx20);
 867             ty               = _mm256_mul_pd(fscal,dy20);
 868             tz               = _mm256_mul_pd(fscal,dz20);
 869
 870             /* Update vectorial force */
 871             fix2             = _mm256_add_pd(fix2,tx);
 872             fiy2             = _mm256_add_pd(fiy2,ty);
 873             fiz2             = _mm256_add_pd(fiz2,tz);
 874
 875             fjx0             = _mm256_add_pd(fjx0,tx);
 876             fjy0             = _mm256_add_pd(fjy0,ty);
 877             fjz0             = _mm256_add_pd(fjz0,tz);
 878
 879             /**************************
 880              * CALCULATE INTERACTIONS *
 881              **************************/
 882
 883             /* COULOMB ELECTROSTATICS */
 884             velec            = _mm256_mul_pd(qq21,rinv21);
 885             felec            = _mm256_mul_pd(velec,rinvsq21);
 886
 887             /* Update potential sum for this i atom from the interaction with this j atom. */
 888             velec            = _mm256_andnot_pd(dummy_mask,velec);
 889             velecsum         = _mm256_add_pd(velecsum,velec);
 890
 891             fscal            = felec;
 892
 893             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
 894
 895             /* Calculate temporary vectorial force */
 896             tx               = _mm256_mul_pd(fscal,dx21);
 897             ty               = _mm256_mul_pd(fscal,dy21);
 898             tz               = _mm256_mul_pd(fscal,dz21);
 899
 900             /* Update vectorial force */
 901             fix2             = _mm256_add_pd(fix2,tx);
 902             fiy2             = _mm256_add_pd(fiy2,ty);
 903             fiz2             = _mm256_add_pd(fiz2,tz);
 904
 905             fjx1             = _mm256_add_pd(fjx1,tx);
 906             fjy1             = _mm256_add_pd(fjy1,ty);
 907             fjz1             = _mm256_add_pd(fjz1,tz);
 908
 909             /**************************
 910              * CALCULATE INTERACTIONS *
 911              **************************/
 912
 913             /* COULOMB ELECTROSTATICS */
 914             velec            = _mm256_mul_pd(qq22,rinv22);
 915             felec            = _mm256_mul_pd(velec,rinvsq22);
 916
 917             /* Update potential sum for this i atom from the interaction with this j atom. */
 918             velec            = _mm256_andnot_pd(dummy_mask,velec);
 919             velecsum         = _mm256_add_pd(velecsum,velec);
 920
 921             fscal            = felec;
 922
 923             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
 924
 925             /* Calculate temporary vectorial force */
 926             tx               = _mm256_mul_pd(fscal,dx22);
 927             ty               = _mm256_mul_pd(fscal,dy22);
 928             tz               = _mm256_mul_pd(fscal,dz22);
 929
 930             /* Update vectorial force */
 931             fix2             = _mm256_add_pd(fix2,tx);
 932             fiy2             = _mm256_add_pd(fiy2,ty);
 933             fiz2             = _mm256_add_pd(fiz2,tz);
 934
 935             fjx2             = _mm256_add_pd(fjx2,tx);
 936             fjy2             = _mm256_add_pd(fjy2,ty);
 937             fjz2             = _mm256_add_pd(fjz2,tz);
 938
 939             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
 940             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
 941             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
 942             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
 943
 944             gmx_mm256_decrement_3rvec_4ptr_swizzle_pd(fjptrA,fjptrB,fjptrC,fjptrD,
 945                                                       fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,fjx2,fjy2,fjz2);
 946
 947             /* Inner loop uses 255 flops */
 948         }
 949
 950         /* End of innermost loop */
 951
 952         gmx_mm256_update_iforce_3atom_swizzle_pd(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,
 953                                                  f+i_coord_offset,fshift+i_shift_offset);
 954
 955         ggid                        = gid[iidx];
 956         /* Update potential energies */
 957         gmx_mm256_update_1pot_pd(velecsum,kernel_data->energygrp_elec+ggid);
 958         gmx_mm256_update_1pot_pd(vvdwsum,kernel_data->energygrp_vdw+ggid);
 959
 960         /* Increment number of inner iterations */
 961         inneriter                  += j_index_end - j_index_start;
 962
 963         /* Outer loop uses 20 flops */
 964     }
 965
 966     /* Increment number of outer iterations */
 967     outeriter        += nri;
 968
 969     /* Update outer/inner flops */
 970
 971     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W3W3_VF,outeriter*20 + inneriter*255);
 972 }
 973 /*
 974  * Gromacs nonbonded kernel:   nb_kernel_ElecCoul_VdwLJ_GeomW3W3_F_avx_256_double
 975  * Electrostatics interaction: Coulomb
 976  * VdW interaction:            LennardJones
 977  * Geometry:                   Water3-Water3
 978  * Calculate force/pot:        Force
 979  */
 980 void
 981 nb_kernel_ElecCoul_VdwLJ_GeomW3W3_F_avx_256_double
 982                     (t_nblist                    * gmx_restrict       nlist,
 983                      rvec                        * gmx_restrict          xx,
 984                      rvec                        * gmx_restrict          ff,
 985                      t_forcerec                  * gmx_restrict          fr,
 986                      t_mdatoms                   * gmx_restrict     mdatoms,
 987                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
 988                      t_nrnb                      * gmx_restrict        nrnb)
 989 {
 990     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
 991      * just 0 for non-waters.
 992      * Suffixes A,B,C,D refer to j loop unrolling done with AVX, e.g. for the four different
 993      * jnr indices corresponding to data put in the four positions in the SIMD register.
 994      */
 995     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
 996     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
 997     int              jnrA,jnrB,jnrC,jnrD;
 998     int              jnrlistA,jnrlistB,jnrlistC,jnrlistD;
 999     int              jnrlistE,jnrlistF,jnrlistG,jnrlistH;
1000     int              j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
1001     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
1002     real             rcutoff_scalar;
1003     real             *shiftvec,*fshift,*x,*f;
1004     real             *fjptrA,*fjptrB,*fjptrC,*fjptrD;
1005     real             scratch[4*DIM];
1006     __m256d          tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
1007     real *           vdwioffsetptr0;
1008     __m256d          ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
1009     real *           vdwioffsetptr1;
1010     __m256d          ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
1011     real *           vdwioffsetptr2;
1012     __m256d          ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
1013     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
1014     __m256d          jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
1015     int              vdwjidx1A,vdwjidx1B,vdwjidx1C,vdwjidx1D;
1016     __m256d          jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
1017     int              vdwjidx2A,vdwjidx2B,vdwjidx2C,vdwjidx2D;
1018     __m256d          jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
1019     __m256d          dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
1020     __m256d          dx01,dy01,dz01,rsq01,rinv01,rinvsq01,r01,qq01,c6_01,c12_01;
1021     __m256d          dx02,dy02,dz02,rsq02,rinv02,rinvsq02,r02,qq02,c6_02,c12_02;
1022     __m256d          dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
1023     __m256d          dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11;
1024     __m256d          dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12;
1025     __m256d          dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
1026     __m256d          dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21;
1027     __m256d          dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22;
1028     __m256d          velec,felec,velecsum,facel,crf,krf,krf2;
1029     real             *charge;
1030     int              nvdwtype;
1031     __m256d          rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
1032     int              *vdwtype;
1033     real             *vdwparam;
1034     __m256d          one_sixth   = _mm256_set1_pd(1.0/6.0);
1035     __m256d          one_twelfth = _mm256_set1_pd(1.0/12.0);
1036     __m256d          dummy_mask,cutoff_mask;
1037     __m128           tmpmask0,tmpmask1;
1038     __m256d          signbit = _mm256_castsi256_pd( _mm256_set1_epi32(0x80000000) );
1039     __m256d          one     = _mm256_set1_pd(1.0);
1040     __m256d          two     = _mm256_set1_pd(2.0);
1041     x                = xx[0];
1042     f                = ff[0];
1043
1044     nri              = nlist->nri;
1045     iinr             = nlist->iinr;
1046     jindex           = nlist->jindex;
1047     jjnr             = nlist->jjnr;
1048     shiftidx         = nlist->shift;
1049     gid              = nlist->gid;
1050     shiftvec         = fr->shift_vec[0];
1051     fshift           = fr->fshift[0];
1052     facel            = _mm256_set1_pd(fr->epsfac);
1053     charge           = mdatoms->chargeA;
1054     nvdwtype         = fr->ntype;
1055     vdwparam         = fr->nbfp;
1056     vdwtype          = mdatoms->typeA;
1057
1058     /* Setup water-specific parameters */
1059     inr              = nlist->iinr[0];
1060     iq0              = _mm256_mul_pd(facel,_mm256_set1_pd(charge[inr+0]));
1061     iq1              = _mm256_mul_pd(facel,_mm256_set1_pd(charge[inr+1]));
1062     iq2              = _mm256_mul_pd(facel,_mm256_set1_pd(charge[inr+2]));
1063     vdwioffsetptr0   = vdwparam+2*nvdwtype*vdwtype[inr+0];
1064
1065     jq0              = _mm256_set1_pd(charge[inr+0]);
1066     jq1              = _mm256_set1_pd(charge[inr+1]);
1067     jq2              = _mm256_set1_pd(charge[inr+2]);
1068     vdwjidx0A        = 2*vdwtype[inr+0];
1069     qq00             = _mm256_mul_pd(iq0,jq0);
1070     c6_00            = _mm256_set1_pd(vdwioffsetptr0[vdwjidx0A]);
1071     c12_00           = _mm256_set1_pd(vdwioffsetptr0[vdwjidx0A+1]);
1072     qq01             = _mm256_mul_pd(iq0,jq1);
1073     qq02             = _mm256_mul_pd(iq0,jq2);
1074     qq10             = _mm256_mul_pd(iq1,jq0);
1075     qq11             = _mm256_mul_pd(iq1,jq1);
1076     qq12             = _mm256_mul_pd(iq1,jq2);
1077     qq20             = _mm256_mul_pd(iq2,jq0);
1078     qq21             = _mm256_mul_pd(iq2,jq1);
1079     qq22             = _mm256_mul_pd(iq2,jq2);
1080
1081     /* Avoid stupid compiler warnings */
1082     jnrA = jnrB = jnrC = jnrD = 0;
1083     j_coord_offsetA = 0;
1084     j_coord_offsetB = 0;
1085     j_coord_offsetC = 0;
1086     j_coord_offsetD = 0;
1087
1088     outeriter        = 0;
1089     inneriter        = 0;
1090
1091     for(iidx=0;iidx<4*DIM;iidx++)
1092     {
1093         scratch[iidx] = 0.0;
1094     }
1095
1096     /* Start outer loop over neighborlists */
1097     for(iidx=0; iidx<nri; iidx++)
1098     {
1099         /* Load shift vector for this list */
1100         i_shift_offset   = DIM*shiftidx[iidx];
1101
1102         /* Load limits for loop over neighbors */
1103         j_index_start    = jindex[iidx];
1104         j_index_end      = jindex[iidx+1];
1105
1106         /* Get outer coordinate index */
1107         inr              = iinr[iidx];
1108         i_coord_offset   = DIM*inr;
1109
1110         /* Load i particle coords and add shift vector */
1111         gmx_mm256_load_shift_and_3rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,
1112                                                     &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2);
1113
1114         fix0             = _mm256_setzero_pd();
1115         fiy0             = _mm256_setzero_pd();
1116         fiz0             = _mm256_setzero_pd();
1117         fix1             = _mm256_setzero_pd();
1118         fiy1             = _mm256_setzero_pd();
1119         fiz1             = _mm256_setzero_pd();
1120         fix2             = _mm256_setzero_pd();
1121         fiy2             = _mm256_setzero_pd();
1122         fiz2             = _mm256_setzero_pd();
1123
1124         /* Start inner kernel loop */
1125         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
1126         {
1127
1128             /* Get j neighbor index, and coordinate index */
1129             jnrA             = jjnr[jidx];
1130             jnrB             = jjnr[jidx+1];
1131             jnrC             = jjnr[jidx+2];
1132             jnrD             = jjnr[jidx+3];
1133             j_coord_offsetA  = DIM*jnrA;
1134             j_coord_offsetB  = DIM*jnrB;
1135             j_coord_offsetC  = DIM*jnrC;
1136             j_coord_offsetD  = DIM*jnrD;
1137
1138             /* load j atom coordinates */
1139             gmx_mm256_load_3rvec_4ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
1140                                                  x+j_coord_offsetC,x+j_coord_offsetD,
1141                                               &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,&jy2,&jz2);
1142
1143             /* Calculate displacement vector */
1144             dx00             = _mm256_sub_pd(ix0,jx0);
1145             dy00             = _mm256_sub_pd(iy0,jy0);
1146             dz00             = _mm256_sub_pd(iz0,jz0);
1147             dx01             = _mm256_sub_pd(ix0,jx1);
1148             dy01             = _mm256_sub_pd(iy0,jy1);
1149             dz01             = _mm256_sub_pd(iz0,jz1);
1150             dx02             = _mm256_sub_pd(ix0,jx2);
1151             dy02             = _mm256_sub_pd(iy0,jy2);
1152             dz02             = _mm256_sub_pd(iz0,jz2);
1153             dx10             = _mm256_sub_pd(ix1,jx0);
1154             dy10             = _mm256_sub_pd(iy1,jy0);
1155             dz10             = _mm256_sub_pd(iz1,jz0);
1156             dx11             = _mm256_sub_pd(ix1,jx1);
1157             dy11             = _mm256_sub_pd(iy1,jy1);
1158             dz11             = _mm256_sub_pd(iz1,jz1);
1159             dx12             = _mm256_sub_pd(ix1,jx2);
1160             dy12             = _mm256_sub_pd(iy1,jy2);
1161             dz12             = _mm256_sub_pd(iz1,jz2);
1162             dx20             = _mm256_sub_pd(ix2,jx0);
1163             dy20             = _mm256_sub_pd(iy2,jy0);
1164             dz20             = _mm256_sub_pd(iz2,jz0);
1165             dx21             = _mm256_sub_pd(ix2,jx1);
1166             dy21             = _mm256_sub_pd(iy2,jy1);
1167             dz21             = _mm256_sub_pd(iz2,jz1);
1168             dx22             = _mm256_sub_pd(ix2,jx2);
1169             dy22             = _mm256_sub_pd(iy2,jy2);
1170             dz22             = _mm256_sub_pd(iz2,jz2);
1171
1172             /* Calculate squared distance and things based on it */
1173             rsq00            = gmx_mm256_calc_rsq_pd(dx00,dy00,dz00);
1174             rsq01            = gmx_mm256_calc_rsq_pd(dx01,dy01,dz01);
1175             rsq02            = gmx_mm256_calc_rsq_pd(dx02,dy02,dz02);
1176             rsq10            = gmx_mm256_calc_rsq_pd(dx10,dy10,dz10);
1177             rsq11            = gmx_mm256_calc_rsq_pd(dx11,dy11,dz11);
1178             rsq12            = gmx_mm256_calc_rsq_pd(dx12,dy12,dz12);
1179             rsq20            = gmx_mm256_calc_rsq_pd(dx20,dy20,dz20);
1180             rsq21            = gmx_mm256_calc_rsq_pd(dx21,dy21,dz21);
1181             rsq22            = gmx_mm256_calc_rsq_pd(dx22,dy22,dz22);
1182
1183             rinv00           = gmx_mm256_invsqrt_pd(rsq00);
1184             rinv01           = gmx_mm256_invsqrt_pd(rsq01);
1185             rinv02           = gmx_mm256_invsqrt_pd(rsq02);
1186             rinv10           = gmx_mm256_invsqrt_pd(rsq10);
1187             rinv11           = gmx_mm256_invsqrt_pd(rsq11);
1188             rinv12           = gmx_mm256_invsqrt_pd(rsq12);
1189             rinv20           = gmx_mm256_invsqrt_pd(rsq20);
1190             rinv21           = gmx_mm256_invsqrt_pd(rsq21);
1191             rinv22           = gmx_mm256_invsqrt_pd(rsq22);
1192
1193             rinvsq00         = _mm256_mul_pd(rinv00,rinv00);
1194             rinvsq01         = _mm256_mul_pd(rinv01,rinv01);
1195             rinvsq02         = _mm256_mul_pd(rinv02,rinv02);
1196             rinvsq10         = _mm256_mul_pd(rinv10,rinv10);
1197             rinvsq11         = _mm256_mul_pd(rinv11,rinv11);
1198             rinvsq12         = _mm256_mul_pd(rinv12,rinv12);
1199             rinvsq20         = _mm256_mul_pd(rinv20,rinv20);
1200             rinvsq21         = _mm256_mul_pd(rinv21,rinv21);
1201             rinvsq22         = _mm256_mul_pd(rinv22,rinv22);
1202
1203             fjx0             = _mm256_setzero_pd();
1204             fjy0             = _mm256_setzero_pd();
1205             fjz0             = _mm256_setzero_pd();
1206             fjx1             = _mm256_setzero_pd();
1207             fjy1             = _mm256_setzero_pd();
1208             fjz1             = _mm256_setzero_pd();
1209             fjx2             = _mm256_setzero_pd();
1210             fjy2             = _mm256_setzero_pd();
1211             fjz2             = _mm256_setzero_pd();
1212
1213             /**************************
1214              * CALCULATE INTERACTIONS *
1215              **************************/
1216
1217             /* COULOMB ELECTROSTATICS */
1218             velec            = _mm256_mul_pd(qq00,rinv00);
1219             felec            = _mm256_mul_pd(velec,rinvsq00);
1220
1221             /* LENNARD-JONES DISPERSION/REPULSION */
1222
1223             rinvsix          = _mm256_mul_pd(_mm256_mul_pd(rinvsq00,rinvsq00),rinvsq00);
1224             fvdw             = _mm256_mul_pd(_mm256_sub_pd(_mm256_mul_pd(c12_00,rinvsix),c6_00),_mm256_mul_pd(rinvsix,rinvsq00));
1225
1226             fscal            = _mm256_add_pd(felec,fvdw);
1227
1228             /* Calculate temporary vectorial force */
1229             tx               = _mm256_mul_pd(fscal,dx00);
1230             ty               = _mm256_mul_pd(fscal,dy00);
1231             tz               = _mm256_mul_pd(fscal,dz00);
1232
1233             /* Update vectorial force */
1234             fix0             = _mm256_add_pd(fix0,tx);
1235             fiy0             = _mm256_add_pd(fiy0,ty);
1236             fiz0             = _mm256_add_pd(fiz0,tz);
1237
1238             fjx0             = _mm256_add_pd(fjx0,tx);
1239             fjy0             = _mm256_add_pd(fjy0,ty);
1240             fjz0             = _mm256_add_pd(fjz0,tz);
1241
1242             /**************************
1243              * CALCULATE INTERACTIONS *
1244              **************************/
1245
1246             /* COULOMB ELECTROSTATICS */
1247             velec            = _mm256_mul_pd(qq01,rinv01);
1248             felec            = _mm256_mul_pd(velec,rinvsq01);
1249
1250             fscal            = felec;
1251
1252             /* Calculate temporary vectorial force */
1253             tx               = _mm256_mul_pd(fscal,dx01);
1254             ty               = _mm256_mul_pd(fscal,dy01);
1255             tz               = _mm256_mul_pd(fscal,dz01);
1256
1257             /* Update vectorial force */
1258             fix0             = _mm256_add_pd(fix0,tx);
1259             fiy0             = _mm256_add_pd(fiy0,ty);
1260             fiz0             = _mm256_add_pd(fiz0,tz);
1261
1262             fjx1             = _mm256_add_pd(fjx1,tx);
1263             fjy1             = _mm256_add_pd(fjy1,ty);
1264             fjz1             = _mm256_add_pd(fjz1,tz);
1265
1266             /**************************
1267              * CALCULATE INTERACTIONS *
1268              **************************/
1269
1270             /* COULOMB ELECTROSTATICS */
1271             velec            = _mm256_mul_pd(qq02,rinv02);
1272             felec            = _mm256_mul_pd(velec,rinvsq02);
1273
1274             fscal            = felec;
1275
1276             /* Calculate temporary vectorial force */
1277             tx               = _mm256_mul_pd(fscal,dx02);
1278             ty               = _mm256_mul_pd(fscal,dy02);
1279             tz               = _mm256_mul_pd(fscal,dz02);
1280
1281             /* Update vectorial force */
1282             fix0             = _mm256_add_pd(fix0,tx);
1283             fiy0             = _mm256_add_pd(fiy0,ty);
1284             fiz0             = _mm256_add_pd(fiz0,tz);
1285
1286             fjx2             = _mm256_add_pd(fjx2,tx);
1287             fjy2             = _mm256_add_pd(fjy2,ty);
1288             fjz2             = _mm256_add_pd(fjz2,tz);
1289
1290             /**************************
1291              * CALCULATE INTERACTIONS *
1292              **************************/
1293
1294             /* COULOMB ELECTROSTATICS */
1295             velec            = _mm256_mul_pd(qq10,rinv10);
1296             felec            = _mm256_mul_pd(velec,rinvsq10);
1297
1298             fscal            = felec;
1299
1300             /* Calculate temporary vectorial force */
1301             tx               = _mm256_mul_pd(fscal,dx10);
1302             ty               = _mm256_mul_pd(fscal,dy10);
1303             tz               = _mm256_mul_pd(fscal,dz10);
1304
1305             /* Update vectorial force */
1306             fix1             = _mm256_add_pd(fix1,tx);
1307             fiy1             = _mm256_add_pd(fiy1,ty);
1308             fiz1             = _mm256_add_pd(fiz1,tz);
1309
1310             fjx0             = _mm256_add_pd(fjx0,tx);
1311             fjy0             = _mm256_add_pd(fjy0,ty);
1312             fjz0             = _mm256_add_pd(fjz0,tz);
1313
1314             /**************************
1315              * CALCULATE INTERACTIONS *
1316              **************************/
1317
1318             /* COULOMB ELECTROSTATICS */
1319             velec            = _mm256_mul_pd(qq11,rinv11);
1320             felec            = _mm256_mul_pd(velec,rinvsq11);
1321
1322             fscal            = felec;
1323
1324             /* Calculate temporary vectorial force */
1325             tx               = _mm256_mul_pd(fscal,dx11);
1326             ty               = _mm256_mul_pd(fscal,dy11);
1327             tz               = _mm256_mul_pd(fscal,dz11);
1328
1329             /* Update vectorial force */
1330             fix1             = _mm256_add_pd(fix1,tx);
1331             fiy1             = _mm256_add_pd(fiy1,ty);
1332             fiz1             = _mm256_add_pd(fiz1,tz);
1333
1334             fjx1             = _mm256_add_pd(fjx1,tx);
1335             fjy1             = _mm256_add_pd(fjy1,ty);
1336             fjz1             = _mm256_add_pd(fjz1,tz);
1337
1338             /**************************
1339              * CALCULATE INTERACTIONS *
1340              **************************/
1341
1342             /* COULOMB ELECTROSTATICS */
1343             velec            = _mm256_mul_pd(qq12,rinv12);
1344             felec            = _mm256_mul_pd(velec,rinvsq12);
1345
1346             fscal            = felec;
1347
1348             /* Calculate temporary vectorial force */
1349             tx               = _mm256_mul_pd(fscal,dx12);
1350             ty               = _mm256_mul_pd(fscal,dy12);
1351             tz               = _mm256_mul_pd(fscal,dz12);
1352
1353             /* Update vectorial force */
1354             fix1             = _mm256_add_pd(fix1,tx);
1355             fiy1             = _mm256_add_pd(fiy1,ty);
1356             fiz1             = _mm256_add_pd(fiz1,tz);
1357
1358             fjx2             = _mm256_add_pd(fjx2,tx);
1359             fjy2             = _mm256_add_pd(fjy2,ty);
1360             fjz2             = _mm256_add_pd(fjz2,tz);
1361
1362             /**************************
1363              * CALCULATE INTERACTIONS *
1364              **************************/
1365
1366             /* COULOMB ELECTROSTATICS */
1367             velec            = _mm256_mul_pd(qq20,rinv20);
1368             felec            = _mm256_mul_pd(velec,rinvsq20);
1369
1370             fscal            = felec;
1371
1372             /* Calculate temporary vectorial force */
1373             tx               = _mm256_mul_pd(fscal,dx20);
1374             ty               = _mm256_mul_pd(fscal,dy20);
1375             tz               = _mm256_mul_pd(fscal,dz20);
1376
1377             /* Update vectorial force */
1378             fix2             = _mm256_add_pd(fix2,tx);
1379             fiy2             = _mm256_add_pd(fiy2,ty);
1380             fiz2             = _mm256_add_pd(fiz2,tz);
1381
1382             fjx0             = _mm256_add_pd(fjx0,tx);
1383             fjy0             = _mm256_add_pd(fjy0,ty);
1384             fjz0             = _mm256_add_pd(fjz0,tz);
1385
1386             /**************************
1387              * CALCULATE INTERACTIONS *
1388              **************************/
1389
1390             /* COULOMB ELECTROSTATICS */
1391             velec            = _mm256_mul_pd(qq21,rinv21);
1392             felec            = _mm256_mul_pd(velec,rinvsq21);
1393
1394             fscal            = felec;
1395
1396             /* Calculate temporary vectorial force */
1397             tx               = _mm256_mul_pd(fscal,dx21);
1398             ty               = _mm256_mul_pd(fscal,dy21);
1399             tz               = _mm256_mul_pd(fscal,dz21);
1400
1401             /* Update vectorial force */
1402             fix2             = _mm256_add_pd(fix2,tx);
1403             fiy2             = _mm256_add_pd(fiy2,ty);
1404             fiz2             = _mm256_add_pd(fiz2,tz);
1405
1406             fjx1             = _mm256_add_pd(fjx1,tx);
1407             fjy1             = _mm256_add_pd(fjy1,ty);
1408             fjz1             = _mm256_add_pd(fjz1,tz);
1409
1410             /**************************
1411              * CALCULATE INTERACTIONS *
1412              **************************/
1413
1414             /* COULOMB ELECTROSTATICS */
1415             velec            = _mm256_mul_pd(qq22,rinv22);
1416             felec            = _mm256_mul_pd(velec,rinvsq22);
1417
1418             fscal            = felec;
1419
1420             /* Calculate temporary vectorial force */
1421             tx               = _mm256_mul_pd(fscal,dx22);
1422             ty               = _mm256_mul_pd(fscal,dy22);
1423             tz               = _mm256_mul_pd(fscal,dz22);
1424
1425             /* Update vectorial force */
1426             fix2             = _mm256_add_pd(fix2,tx);
1427             fiy2             = _mm256_add_pd(fiy2,ty);
1428             fiz2             = _mm256_add_pd(fiz2,tz);
1429
1430             fjx2             = _mm256_add_pd(fjx2,tx);
1431             fjy2             = _mm256_add_pd(fjy2,ty);
1432             fjz2             = _mm256_add_pd(fjz2,tz);
1433
1434             fjptrA             = f+j_coord_offsetA;
1435             fjptrB             = f+j_coord_offsetB;
1436             fjptrC             = f+j_coord_offsetC;
1437             fjptrD             = f+j_coord_offsetD;
1438
1439             gmx_mm256_decrement_3rvec_4ptr_swizzle_pd(fjptrA,fjptrB,fjptrC,fjptrD,
1440                                                       fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,fjx2,fjy2,fjz2);
1441
1442             /* Inner loop uses 241 flops */
1443         }
1444
1445         if(jidx<j_index_end)
1446         {
1447
1448             /* Get j neighbor index, and coordinate index */
1449             jnrlistA         = jjnr[jidx];
1450             jnrlistB         = jjnr[jidx+1];
1451             jnrlistC         = jjnr[jidx+2];
1452             jnrlistD         = jjnr[jidx+3];
1453             /* Sign of each element will be negative for non-real atoms.
1454              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
1455              * so use it as val = _mm_andnot_pd(mask,val) to clear dummy entries.
1456              */
1457             tmpmask0 = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
1458
1459             tmpmask1 = _mm_permute_ps(tmpmask0,_GMX_MM_PERMUTE(3,3,2,2));
1460             tmpmask0 = _mm_permute_ps(tmpmask0,_GMX_MM_PERMUTE(1,1,0,0));
1461             dummy_mask = _mm256_castps_pd(gmx_mm256_set_m128(tmpmask1,tmpmask0));
1462
1463             jnrA       = (jnrlistA>=0) ? jnrlistA : 0;
1464             jnrB       = (jnrlistB>=0) ? jnrlistB : 0;
1465             jnrC       = (jnrlistC>=0) ? jnrlistC : 0;
1466             jnrD       = (jnrlistD>=0) ? jnrlistD : 0;
1467             j_coord_offsetA  = DIM*jnrA;
1468             j_coord_offsetB  = DIM*jnrB;
1469             j_coord_offsetC  = DIM*jnrC;
1470             j_coord_offsetD  = DIM*jnrD;
1471
1472             /* load j atom coordinates */
1473             gmx_mm256_load_3rvec_4ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
1474                                                  x+j_coord_offsetC,x+j_coord_offsetD,
1475                                               &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,&jy2,&jz2);
1476
1477             /* Calculate displacement vector */
1478             dx00             = _mm256_sub_pd(ix0,jx0);
1479             dy00             = _mm256_sub_pd(iy0,jy0);
1480             dz00             = _mm256_sub_pd(iz0,jz0);
1481             dx01             = _mm256_sub_pd(ix0,jx1);
1482             dy01             = _mm256_sub_pd(iy0,jy1);
1483             dz01             = _mm256_sub_pd(iz0,jz1);
1484             dx02             = _mm256_sub_pd(ix0,jx2);
1485             dy02             = _mm256_sub_pd(iy0,jy2);
1486             dz02             = _mm256_sub_pd(iz0,jz2);
1487             dx10             = _mm256_sub_pd(ix1,jx0);
1488             dy10             = _mm256_sub_pd(iy1,jy0);
1489             dz10             = _mm256_sub_pd(iz1,jz0);
1490             dx11             = _mm256_sub_pd(ix1,jx1);
1491             dy11             = _mm256_sub_pd(iy1,jy1);
1492             dz11             = _mm256_sub_pd(iz1,jz1);
1493             dx12             = _mm256_sub_pd(ix1,jx2);
1494             dy12             = _mm256_sub_pd(iy1,jy2);
1495             dz12             = _mm256_sub_pd(iz1,jz2);
1496             dx20             = _mm256_sub_pd(ix2,jx0);
1497             dy20             = _mm256_sub_pd(iy2,jy0);
1498             dz20             = _mm256_sub_pd(iz2,jz0);
1499             dx21             = _mm256_sub_pd(ix2,jx1);
1500             dy21             = _mm256_sub_pd(iy2,jy1);
1501             dz21             = _mm256_sub_pd(iz2,jz1);
1502             dx22             = _mm256_sub_pd(ix2,jx2);
1503             dy22             = _mm256_sub_pd(iy2,jy2);
1504             dz22             = _mm256_sub_pd(iz2,jz2);
1505
1506             /* Calculate squared distance and things based on it */
1507             rsq00            = gmx_mm256_calc_rsq_pd(dx00,dy00,dz00);
1508             rsq01            = gmx_mm256_calc_rsq_pd(dx01,dy01,dz01);
1509             rsq02            = gmx_mm256_calc_rsq_pd(dx02,dy02,dz02);
1510             rsq10            = gmx_mm256_calc_rsq_pd(dx10,dy10,dz10);
1511             rsq11            = gmx_mm256_calc_rsq_pd(dx11,dy11,dz11);
1512             rsq12            = gmx_mm256_calc_rsq_pd(dx12,dy12,dz12);
1513             rsq20            = gmx_mm256_calc_rsq_pd(dx20,dy20,dz20);
1514             rsq21            = gmx_mm256_calc_rsq_pd(dx21,dy21,dz21);
1515             rsq22            = gmx_mm256_calc_rsq_pd(dx22,dy22,dz22);
1516
1517             rinv00           = gmx_mm256_invsqrt_pd(rsq00);
1518             rinv01           = gmx_mm256_invsqrt_pd(rsq01);
1519             rinv02           = gmx_mm256_invsqrt_pd(rsq02);
1520             rinv10           = gmx_mm256_invsqrt_pd(rsq10);
1521             rinv11           = gmx_mm256_invsqrt_pd(rsq11);
1522             rinv12           = gmx_mm256_invsqrt_pd(rsq12);
1523             rinv20           = gmx_mm256_invsqrt_pd(rsq20);
1524             rinv21           = gmx_mm256_invsqrt_pd(rsq21);
1525             rinv22           = gmx_mm256_invsqrt_pd(rsq22);
1526
1527             rinvsq00         = _mm256_mul_pd(rinv00,rinv00);
1528             rinvsq01         = _mm256_mul_pd(rinv01,rinv01);
1529             rinvsq02         = _mm256_mul_pd(rinv02,rinv02);
1530             rinvsq10         = _mm256_mul_pd(rinv10,rinv10);
1531             rinvsq11         = _mm256_mul_pd(rinv11,rinv11);
1532             rinvsq12         = _mm256_mul_pd(rinv12,rinv12);
1533             rinvsq20         = _mm256_mul_pd(rinv20,rinv20);
1534             rinvsq21         = _mm256_mul_pd(rinv21,rinv21);
1535             rinvsq22         = _mm256_mul_pd(rinv22,rinv22);
1536
1537             fjx0             = _mm256_setzero_pd();
1538             fjy0             = _mm256_setzero_pd();
1539             fjz0             = _mm256_setzero_pd();
1540             fjx1             = _mm256_setzero_pd();
1541             fjy1             = _mm256_setzero_pd();
1542             fjz1             = _mm256_setzero_pd();
1543             fjx2             = _mm256_setzero_pd();
1544             fjy2             = _mm256_setzero_pd();
1545             fjz2             = _mm256_setzero_pd();
1546
1547             /**************************
1548              * CALCULATE INTERACTIONS *
1549              **************************/
1550
1551             /* COULOMB ELECTROSTATICS */
1552             velec            = _mm256_mul_pd(qq00,rinv00);
1553             felec            = _mm256_mul_pd(velec,rinvsq00);
1554
1555             /* LENNARD-JONES DISPERSION/REPULSION */
1556
1557             rinvsix          = _mm256_mul_pd(_mm256_mul_pd(rinvsq00,rinvsq00),rinvsq00);
1558             fvdw             = _mm256_mul_pd(_mm256_sub_pd(_mm256_mul_pd(c12_00,rinvsix),c6_00),_mm256_mul_pd(rinvsix,rinvsq00));
1559
1560             fscal            = _mm256_add_pd(felec,fvdw);
1561
1562             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
1563
1564             /* Calculate temporary vectorial force */
1565             tx               = _mm256_mul_pd(fscal,dx00);
1566             ty               = _mm256_mul_pd(fscal,dy00);
1567             tz               = _mm256_mul_pd(fscal,dz00);
1568
1569             /* Update vectorial force */
1570             fix0             = _mm256_add_pd(fix0,tx);
1571             fiy0             = _mm256_add_pd(fiy0,ty);
1572             fiz0             = _mm256_add_pd(fiz0,tz);
1573
1574             fjx0             = _mm256_add_pd(fjx0,tx);
1575             fjy0             = _mm256_add_pd(fjy0,ty);
1576             fjz0             = _mm256_add_pd(fjz0,tz);
1577
1578             /**************************
1579              * CALCULATE INTERACTIONS *
1580              **************************/
1581
1582             /* COULOMB ELECTROSTATICS */
1583             velec            = _mm256_mul_pd(qq01,rinv01);
1584             felec            = _mm256_mul_pd(velec,rinvsq01);
1585
1586             fscal            = felec;
1587
1588             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
1589
1590             /* Calculate temporary vectorial force */
1591             tx               = _mm256_mul_pd(fscal,dx01);
1592             ty               = _mm256_mul_pd(fscal,dy01);
1593             tz               = _mm256_mul_pd(fscal,dz01);
1594
1595             /* Update vectorial force */
1596             fix0             = _mm256_add_pd(fix0,tx);
1597             fiy0             = _mm256_add_pd(fiy0,ty);
1598             fiz0             = _mm256_add_pd(fiz0,tz);
1599
1600             fjx1             = _mm256_add_pd(fjx1,tx);
1601             fjy1             = _mm256_add_pd(fjy1,ty);
1602             fjz1             = _mm256_add_pd(fjz1,tz);
1603
1604             /**************************
1605              * CALCULATE INTERACTIONS *
1606              **************************/
1607
1608             /* COULOMB ELECTROSTATICS */
1609             velec            = _mm256_mul_pd(qq02,rinv02);
1610             felec            = _mm256_mul_pd(velec,rinvsq02);
1611
1612             fscal            = felec;
1613
1614             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
1615
1616             /* Calculate temporary vectorial force */
1617             tx               = _mm256_mul_pd(fscal,dx02);
1618             ty               = _mm256_mul_pd(fscal,dy02);
1619             tz               = _mm256_mul_pd(fscal,dz02);
1620
1621             /* Update vectorial force */
1622             fix0             = _mm256_add_pd(fix0,tx);
1623             fiy0             = _mm256_add_pd(fiy0,ty);
1624             fiz0             = _mm256_add_pd(fiz0,tz);
1625
1626             fjx2             = _mm256_add_pd(fjx2,tx);
1627             fjy2             = _mm256_add_pd(fjy2,ty);
1628             fjz2             = _mm256_add_pd(fjz2,tz);
1629
1630             /**************************
1631              * CALCULATE INTERACTIONS *
1632              **************************/
1633
1634             /* COULOMB ELECTROSTATICS */
1635             velec            = _mm256_mul_pd(qq10,rinv10);
1636             felec            = _mm256_mul_pd(velec,rinvsq10);
1637
1638             fscal            = felec;
1639
1640             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
1641
1642             /* Calculate temporary vectorial force */
1643             tx               = _mm256_mul_pd(fscal,dx10);
1644             ty               = _mm256_mul_pd(fscal,dy10);
1645             tz               = _mm256_mul_pd(fscal,dz10);
1646
1647             /* Update vectorial force */
1648             fix1             = _mm256_add_pd(fix1,tx);
1649             fiy1             = _mm256_add_pd(fiy1,ty);
1650             fiz1             = _mm256_add_pd(fiz1,tz);
1651
1652             fjx0             = _mm256_add_pd(fjx0,tx);
1653             fjy0             = _mm256_add_pd(fjy0,ty);
1654             fjz0             = _mm256_add_pd(fjz0,tz);
1655
1656             /**************************
1657              * CALCULATE INTERACTIONS *
1658              **************************/
1659
1660             /* COULOMB ELECTROSTATICS */
1661             velec            = _mm256_mul_pd(qq11,rinv11);
1662             felec            = _mm256_mul_pd(velec,rinvsq11);
1663
1664             fscal            = felec;
1665
1666             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
1667
1668             /* Calculate temporary vectorial force */
1669             tx               = _mm256_mul_pd(fscal,dx11);
1670             ty               = _mm256_mul_pd(fscal,dy11);
1671             tz               = _mm256_mul_pd(fscal,dz11);
1672
1673             /* Update vectorial force */
1674             fix1             = _mm256_add_pd(fix1,tx);
1675             fiy1             = _mm256_add_pd(fiy1,ty);
1676             fiz1             = _mm256_add_pd(fiz1,tz);
1677
1678             fjx1             = _mm256_add_pd(fjx1,tx);
1679             fjy1             = _mm256_add_pd(fjy1,ty);
1680             fjz1             = _mm256_add_pd(fjz1,tz);
1681
1682             /**************************
1683              * CALCULATE INTERACTIONS *
1684              **************************/
1685
1686             /* COULOMB ELECTROSTATICS */
1687             velec            = _mm256_mul_pd(qq12,rinv12);
1688             felec            = _mm256_mul_pd(velec,rinvsq12);
1689
1690             fscal            = felec;
1691
1692             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
1693
1694             /* Calculate temporary vectorial force */
1695             tx               = _mm256_mul_pd(fscal,dx12);
1696             ty               = _mm256_mul_pd(fscal,dy12);
1697             tz               = _mm256_mul_pd(fscal,dz12);
1698
1699             /* Update vectorial force */
1700             fix1             = _mm256_add_pd(fix1,tx);
1701             fiy1             = _mm256_add_pd(fiy1,ty);
1702             fiz1             = _mm256_add_pd(fiz1,tz);
1703
1704             fjx2             = _mm256_add_pd(fjx2,tx);
1705             fjy2             = _mm256_add_pd(fjy2,ty);
1706             fjz2             = _mm256_add_pd(fjz2,tz);
1707
1708             /**************************
1709              * CALCULATE INTERACTIONS *
1710              **************************/
1711
1712             /* COULOMB ELECTROSTATICS */
1713             velec            = _mm256_mul_pd(qq20,rinv20);
1714             felec            = _mm256_mul_pd(velec,rinvsq20);
1715
1716             fscal            = felec;
1717
1718             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
1719
1720             /* Calculate temporary vectorial force */
1721             tx               = _mm256_mul_pd(fscal,dx20);
1722             ty               = _mm256_mul_pd(fscal,dy20);
1723             tz               = _mm256_mul_pd(fscal,dz20);
1724
1725             /* Update vectorial force */
1726             fix2             = _mm256_add_pd(fix2,tx);
1727             fiy2             = _mm256_add_pd(fiy2,ty);
1728             fiz2             = _mm256_add_pd(fiz2,tz);
1729
1730             fjx0             = _mm256_add_pd(fjx0,tx);
1731             fjy0             = _mm256_add_pd(fjy0,ty);
1732             fjz0             = _mm256_add_pd(fjz0,tz);
1733
1734             /**************************
1735              * CALCULATE INTERACTIONS *
1736              **************************/
1737
1738             /* COULOMB ELECTROSTATICS */
1739             velec            = _mm256_mul_pd(qq21,rinv21);
1740             felec            = _mm256_mul_pd(velec,rinvsq21);
1741
1742             fscal            = felec;
1743
1744             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
1745
1746             /* Calculate temporary vectorial force */
1747             tx               = _mm256_mul_pd(fscal,dx21);
1748             ty               = _mm256_mul_pd(fscal,dy21);
1749             tz               = _mm256_mul_pd(fscal,dz21);
1750
1751             /* Update vectorial force */
1752             fix2             = _mm256_add_pd(fix2,tx);
1753             fiy2             = _mm256_add_pd(fiy2,ty);
1754             fiz2             = _mm256_add_pd(fiz2,tz);
1755
1756             fjx1             = _mm256_add_pd(fjx1,tx);
1757             fjy1             = _mm256_add_pd(fjy1,ty);
1758             fjz1             = _mm256_add_pd(fjz1,tz);
1759
1760             /**************************
1761              * CALCULATE INTERACTIONS *
1762              **************************/
1763
1764             /* COULOMB ELECTROSTATICS */
1765             velec            = _mm256_mul_pd(qq22,rinv22);
1766             felec            = _mm256_mul_pd(velec,rinvsq22);
1767
1768             fscal            = felec;
1769
1770             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
1771
1772             /* Calculate temporary vectorial force */
1773             tx               = _mm256_mul_pd(fscal,dx22);
1774             ty               = _mm256_mul_pd(fscal,dy22);
1775             tz               = _mm256_mul_pd(fscal,dz22);
1776
1777             /* Update vectorial force */
1778             fix2             = _mm256_add_pd(fix2,tx);
1779             fiy2             = _mm256_add_pd(fiy2,ty);
1780             fiz2             = _mm256_add_pd(fiz2,tz);
1781
1782             fjx2             = _mm256_add_pd(fjx2,tx);
1783             fjy2             = _mm256_add_pd(fjy2,ty);
1784             fjz2             = _mm256_add_pd(fjz2,tz);
1785
1786             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
1787             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
1788             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
1789             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
1790
1791             gmx_mm256_decrement_3rvec_4ptr_swizzle_pd(fjptrA,fjptrB,fjptrC,fjptrD,
1792                                                       fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,fjx2,fjy2,fjz2);
1793
1794             /* Inner loop uses 241 flops */
1795         }
1796
1797         /* End of innermost loop */
1798
1799         gmx_mm256_update_iforce_3atom_swizzle_pd(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,
1800                                                  f+i_coord_offset,fshift+i_shift_offset);
1801
1802         /* Increment number of inner iterations */
1803         inneriter                  += j_index_end - j_index_start;
1804
1805         /* Outer loop uses 18 flops */
1806     }
1807
1808     /* Increment number of outer iterations */
1809     outeriter        += nri;
1810
1811     /* Update outer/inner flops */
1812
1813     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W3W3_F,outeriter*18 + inneriter*241);
1814 }