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