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

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*
  36  * Note: this file was generated by the GROMACS avx_256_double kernel generator.
  37  */
  38 #include "config.h"
  39
  40 #include <math.h>
  41
  42 #include "../nb_kernel.h"
  43 #include "gromacs/legacyheaders/types/simple.h"
  44 #include "gromacs/math/vec.h"
  45 #include "gromacs/legacyheaders/nrnb.h"
  46
  47 #include "gromacs/simd/math_x86_avx_256_double.h"
  48 #include "kernelutil_x86_avx_256_double.h"
  49
  50 /*
  51  * Gromacs nonbonded kernel:   nb_kernel_ElecEwSh_VdwNone_GeomW4P1_VF_avx_256_double
  52  * Electrostatics interaction: Ewald
  53  * VdW interaction:            None
  54  * Geometry:                   Water4-Particle
  55  * Calculate force/pot:        PotentialAndForce
  56  */
  57 void
  58 nb_kernel_ElecEwSh_VdwNone_GeomW4P1_VF_avx_256_double
  59                     (t_nblist                    * gmx_restrict       nlist,
  60                      rvec                        * gmx_restrict          xx,
  61                      rvec                        * gmx_restrict          ff,
  62                      t_forcerec                  * gmx_restrict          fr,
  63                      t_mdatoms                   * gmx_restrict     mdatoms,
  64                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
  65                      t_nrnb                      * gmx_restrict        nrnb)
  66 {
  67     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
  68      * just 0 for non-waters.
  69      * Suffixes A,B,C,D refer to j loop unrolling done with AVX, e.g. for the four different
  70      * jnr indices corresponding to data put in the four positions in the SIMD register.
  71      */
  72     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
  73     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
  74     int              jnrA,jnrB,jnrC,jnrD;
  75     int              jnrlistA,jnrlistB,jnrlistC,jnrlistD;
  76     int              jnrlistE,jnrlistF,jnrlistG,jnrlistH;
  77     int              j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
  78     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
  79     real             rcutoff_scalar;
  80     real             *shiftvec,*fshift,*x,*f;
  81     real             *fjptrA,*fjptrB,*fjptrC,*fjptrD;
  82     real             scratch[4*DIM];
  83     __m256d          tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
  84     real *           vdwioffsetptr1;
  85     __m256d          ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
  86     real *           vdwioffsetptr2;
  87     __m256d          ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
  88     real *           vdwioffsetptr3;
  89     __m256d          ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
  90     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
  91     __m256d          jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
  92     __m256d          dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
  93     __m256d          dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
  94     __m256d          dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
  95     __m256d          velec,felec,velecsum,facel,crf,krf,krf2;
  96     real             *charge;
  97     __m128i          ewitab;
  98     __m256d          ewtabscale,eweps,sh_ewald,ewrt,ewtabhalfspace,ewtabF,ewtabFn,ewtabD,ewtabV;
  99     __m256d          beta,beta2,beta3,zeta2,pmecorrF,pmecorrV,rinv3;
 100     real             *ewtab;
 101     __m256d          dummy_mask,cutoff_mask;
 102     __m128           tmpmask0,tmpmask1;
 103     __m256d          signbit = _mm256_castsi256_pd( _mm256_set1_epi32(0x80000000) );
 104     __m256d          one     = _mm256_set1_pd(1.0);
 105     __m256d          two     = _mm256_set1_pd(2.0);
 106     x                = xx[0];
 107     f                = ff[0];
 108
 109     nri              = nlist->nri;
 110     iinr             = nlist->iinr;
 111     jindex           = nlist->jindex;
 112     jjnr             = nlist->jjnr;
 113     shiftidx         = nlist->shift;
 114     gid              = nlist->gid;
 115     shiftvec         = fr->shift_vec[0];
 116     fshift           = fr->fshift[0];
 117     facel            = _mm256_set1_pd(fr->epsfac);
 118     charge           = mdatoms->chargeA;
 119
 120     sh_ewald         = _mm256_set1_pd(fr->ic->sh_ewald);
 121     beta             = _mm256_set1_pd(fr->ic->ewaldcoeff_q);
 122     beta2            = _mm256_mul_pd(beta,beta);
 123     beta3            = _mm256_mul_pd(beta,beta2);
 124
 125     ewtab            = fr->ic->tabq_coul_FDV0;
 126     ewtabscale       = _mm256_set1_pd(fr->ic->tabq_scale);
 127     ewtabhalfspace   = _mm256_set1_pd(0.5/fr->ic->tabq_scale);
 128
 129     /* Setup water-specific parameters */
 130     inr              = nlist->iinr[0];
 131     iq1              = _mm256_mul_pd(facel,_mm256_set1_pd(charge[inr+1]));
 132     iq2              = _mm256_mul_pd(facel,_mm256_set1_pd(charge[inr+2]));
 133     iq3              = _mm256_mul_pd(facel,_mm256_set1_pd(charge[inr+3]));
 134
 135     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
 136     rcutoff_scalar   = fr->rcoulomb;
 137     rcutoff          = _mm256_set1_pd(rcutoff_scalar);
 138     rcutoff2         = _mm256_mul_pd(rcutoff,rcutoff);
 139
 140     /* Avoid stupid compiler warnings */
 141     jnrA = jnrB = jnrC = jnrD = 0;
 142     j_coord_offsetA = 0;
 143     j_coord_offsetB = 0;
 144     j_coord_offsetC = 0;
 145     j_coord_offsetD = 0;
 146
 147     outeriter        = 0;
 148     inneriter        = 0;
 149
 150     for(iidx=0;iidx<4*DIM;iidx++)
 151     {
 152         scratch[iidx] = 0.0;
 153     }
 154
 155     /* Start outer loop over neighborlists */
 156     for(iidx=0; iidx<nri; iidx++)
 157     {
 158         /* Load shift vector for this list */
 159         i_shift_offset   = DIM*shiftidx[iidx];
 160
 161         /* Load limits for loop over neighbors */
 162         j_index_start    = jindex[iidx];
 163         j_index_end      = jindex[iidx+1];
 164
 165         /* Get outer coordinate index */
 166         inr              = iinr[iidx];
 167         i_coord_offset   = DIM*inr;
 168
 169         /* Load i particle coords and add shift vector */
 170         gmx_mm256_load_shift_and_3rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset+DIM,
 171                                                     &ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
 172
 173         fix1             = _mm256_setzero_pd();
 174         fiy1             = _mm256_setzero_pd();
 175         fiz1             = _mm256_setzero_pd();
 176         fix2             = _mm256_setzero_pd();
 177         fiy2             = _mm256_setzero_pd();
 178         fiz2             = _mm256_setzero_pd();
 179         fix3             = _mm256_setzero_pd();
 180         fiy3             = _mm256_setzero_pd();
 181         fiz3             = _mm256_setzero_pd();
 182
 183         /* Reset potential sums */
 184         velecsum         = _mm256_setzero_pd();
 185
 186         /* Start inner kernel loop */
 187         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
 188         {
 189
 190             /* Get j neighbor index, and coordinate index */
 191             jnrA             = jjnr[jidx];
 192             jnrB             = jjnr[jidx+1];
 193             jnrC             = jjnr[jidx+2];
 194             jnrD             = jjnr[jidx+3];
 195             j_coord_offsetA  = DIM*jnrA;
 196             j_coord_offsetB  = DIM*jnrB;
 197             j_coord_offsetC  = DIM*jnrC;
 198             j_coord_offsetD  = DIM*jnrD;
 199
 200             /* load j atom coordinates */
 201             gmx_mm256_load_1rvec_4ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
 202                                                  x+j_coord_offsetC,x+j_coord_offsetD,
 203                                                  &jx0,&jy0,&jz0);
 204
 205             /* Calculate displacement vector */
 206             dx10             = _mm256_sub_pd(ix1,jx0);
 207             dy10             = _mm256_sub_pd(iy1,jy0);
 208             dz10             = _mm256_sub_pd(iz1,jz0);
 209             dx20             = _mm256_sub_pd(ix2,jx0);
 210             dy20             = _mm256_sub_pd(iy2,jy0);
 211             dz20             = _mm256_sub_pd(iz2,jz0);
 212             dx30             = _mm256_sub_pd(ix3,jx0);
 213             dy30             = _mm256_sub_pd(iy3,jy0);
 214             dz30             = _mm256_sub_pd(iz3,jz0);
 215
 216             /* Calculate squared distance and things based on it */
 217             rsq10            = gmx_mm256_calc_rsq_pd(dx10,dy10,dz10);
 218             rsq20            = gmx_mm256_calc_rsq_pd(dx20,dy20,dz20);
 219             rsq30            = gmx_mm256_calc_rsq_pd(dx30,dy30,dz30);
 220
 221             rinv10           = gmx_mm256_invsqrt_pd(rsq10);
 222             rinv20           = gmx_mm256_invsqrt_pd(rsq20);
 223             rinv30           = gmx_mm256_invsqrt_pd(rsq30);
 224
 225             rinvsq10         = _mm256_mul_pd(rinv10,rinv10);
 226             rinvsq20         = _mm256_mul_pd(rinv20,rinv20);
 227             rinvsq30         = _mm256_mul_pd(rinv30,rinv30);
 228
 229             /* Load parameters for j particles */
 230             jq0              = gmx_mm256_load_4real_swizzle_pd(charge+jnrA+0,charge+jnrB+0,
 231                                                                  charge+jnrC+0,charge+jnrD+0);
 232
 233             fjx0             = _mm256_setzero_pd();
 234             fjy0             = _mm256_setzero_pd();
 235             fjz0             = _mm256_setzero_pd();
 236
 237             /**************************
 238              * CALCULATE INTERACTIONS *
 239              **************************/
 240
 241             if (gmx_mm256_any_lt(rsq10,rcutoff2))
 242             {
 243
 244             r10              = _mm256_mul_pd(rsq10,rinv10);
 245
 246             /* Compute parameters for interactions between i and j atoms */
 247             qq10             = _mm256_mul_pd(iq1,jq0);
 248
 249             /* EWALD ELECTROSTATICS */
 250
 251             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 252             ewrt             = _mm256_mul_pd(r10,ewtabscale);
 253             ewitab           = _mm256_cvttpd_epi32(ewrt);
 254             eweps            = _mm256_sub_pd(ewrt,_mm256_round_pd(ewrt, _MM_FROUND_FLOOR));
 255             ewitab           = _mm_slli_epi32(ewitab,2);
 256             ewtabF           = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 257             ewtabD           = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 258             ewtabV           = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,2) );
 259             ewtabFn          = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,3) );
 260             GMX_MM256_FULLTRANSPOSE4_PD(ewtabF,ewtabD,ewtabV,ewtabFn);
 261             felec            = _mm256_add_pd(ewtabF,_mm256_mul_pd(eweps,ewtabD));
 262             velec            = _mm256_sub_pd(ewtabV,_mm256_mul_pd(_mm256_mul_pd(ewtabhalfspace,eweps),_mm256_add_pd(ewtabF,felec)));
 263             velec            = _mm256_mul_pd(qq10,_mm256_sub_pd(_mm256_sub_pd(rinv10,sh_ewald),velec));
 264             felec            = _mm256_mul_pd(_mm256_mul_pd(qq10,rinv10),_mm256_sub_pd(rinvsq10,felec));
 265
 266             cutoff_mask      = _mm256_cmp_pd(rsq10,rcutoff2,_CMP_LT_OQ);
 267
 268             /* Update potential sum for this i atom from the interaction with this j atom. */
 269             velec            = _mm256_and_pd(velec,cutoff_mask);
 270             velecsum         = _mm256_add_pd(velecsum,velec);
 271
 272             fscal            = felec;
 273
 274             fscal            = _mm256_and_pd(fscal,cutoff_mask);
 275
 276             /* Calculate temporary vectorial force */
 277             tx               = _mm256_mul_pd(fscal,dx10);
 278             ty               = _mm256_mul_pd(fscal,dy10);
 279             tz               = _mm256_mul_pd(fscal,dz10);
 280
 281             /* Update vectorial force */
 282             fix1             = _mm256_add_pd(fix1,tx);
 283             fiy1             = _mm256_add_pd(fiy1,ty);
 284             fiz1             = _mm256_add_pd(fiz1,tz);
 285
 286             fjx0             = _mm256_add_pd(fjx0,tx);
 287             fjy0             = _mm256_add_pd(fjy0,ty);
 288             fjz0             = _mm256_add_pd(fjz0,tz);
 289
 290             }
 291
 292             /**************************
 293              * CALCULATE INTERACTIONS *
 294              **************************/
 295
 296             if (gmx_mm256_any_lt(rsq20,rcutoff2))
 297             {
 298
 299             r20              = _mm256_mul_pd(rsq20,rinv20);
 300
 301             /* Compute parameters for interactions between i and j atoms */
 302             qq20             = _mm256_mul_pd(iq2,jq0);
 303
 304             /* EWALD ELECTROSTATICS */
 305
 306             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 307             ewrt             = _mm256_mul_pd(r20,ewtabscale);
 308             ewitab           = _mm256_cvttpd_epi32(ewrt);
 309             eweps            = _mm256_sub_pd(ewrt,_mm256_round_pd(ewrt, _MM_FROUND_FLOOR));
 310             ewitab           = _mm_slli_epi32(ewitab,2);
 311             ewtabF           = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 312             ewtabD           = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 313             ewtabV           = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,2) );
 314             ewtabFn          = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,3) );
 315             GMX_MM256_FULLTRANSPOSE4_PD(ewtabF,ewtabD,ewtabV,ewtabFn);
 316             felec            = _mm256_add_pd(ewtabF,_mm256_mul_pd(eweps,ewtabD));
 317             velec            = _mm256_sub_pd(ewtabV,_mm256_mul_pd(_mm256_mul_pd(ewtabhalfspace,eweps),_mm256_add_pd(ewtabF,felec)));
 318             velec            = _mm256_mul_pd(qq20,_mm256_sub_pd(_mm256_sub_pd(rinv20,sh_ewald),velec));
 319             felec            = _mm256_mul_pd(_mm256_mul_pd(qq20,rinv20),_mm256_sub_pd(rinvsq20,felec));
 320
 321             cutoff_mask      = _mm256_cmp_pd(rsq20,rcutoff2,_CMP_LT_OQ);
 322
 323             /* Update potential sum for this i atom from the interaction with this j atom. */
 324             velec            = _mm256_and_pd(velec,cutoff_mask);
 325             velecsum         = _mm256_add_pd(velecsum,velec);
 326
 327             fscal            = felec;
 328
 329             fscal            = _mm256_and_pd(fscal,cutoff_mask);
 330
 331             /* Calculate temporary vectorial force */
 332             tx               = _mm256_mul_pd(fscal,dx20);
 333             ty               = _mm256_mul_pd(fscal,dy20);
 334             tz               = _mm256_mul_pd(fscal,dz20);
 335
 336             /* Update vectorial force */
 337             fix2             = _mm256_add_pd(fix2,tx);
 338             fiy2             = _mm256_add_pd(fiy2,ty);
 339             fiz2             = _mm256_add_pd(fiz2,tz);
 340
 341             fjx0             = _mm256_add_pd(fjx0,tx);
 342             fjy0             = _mm256_add_pd(fjy0,ty);
 343             fjz0             = _mm256_add_pd(fjz0,tz);
 344
 345             }
 346
 347             /**************************
 348              * CALCULATE INTERACTIONS *
 349              **************************/
 350
 351             if (gmx_mm256_any_lt(rsq30,rcutoff2))
 352             {
 353
 354             r30              = _mm256_mul_pd(rsq30,rinv30);
 355
 356             /* Compute parameters for interactions between i and j atoms */
 357             qq30             = _mm256_mul_pd(iq3,jq0);
 358
 359             /* EWALD ELECTROSTATICS */
 360
 361             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 362             ewrt             = _mm256_mul_pd(r30,ewtabscale);
 363             ewitab           = _mm256_cvttpd_epi32(ewrt);
 364             eweps            = _mm256_sub_pd(ewrt,_mm256_round_pd(ewrt, _MM_FROUND_FLOOR));
 365             ewitab           = _mm_slli_epi32(ewitab,2);
 366             ewtabF           = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 367             ewtabD           = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 368             ewtabV           = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,2) );
 369             ewtabFn          = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,3) );
 370             GMX_MM256_FULLTRANSPOSE4_PD(ewtabF,ewtabD,ewtabV,ewtabFn);
 371             felec            = _mm256_add_pd(ewtabF,_mm256_mul_pd(eweps,ewtabD));
 372             velec            = _mm256_sub_pd(ewtabV,_mm256_mul_pd(_mm256_mul_pd(ewtabhalfspace,eweps),_mm256_add_pd(ewtabF,felec)));
 373             velec            = _mm256_mul_pd(qq30,_mm256_sub_pd(_mm256_sub_pd(rinv30,sh_ewald),velec));
 374             felec            = _mm256_mul_pd(_mm256_mul_pd(qq30,rinv30),_mm256_sub_pd(rinvsq30,felec));
 375
 376             cutoff_mask      = _mm256_cmp_pd(rsq30,rcutoff2,_CMP_LT_OQ);
 377
 378             /* Update potential sum for this i atom from the interaction with this j atom. */
 379             velec            = _mm256_and_pd(velec,cutoff_mask);
 380             velecsum         = _mm256_add_pd(velecsum,velec);
 381
 382             fscal            = felec;
 383
 384             fscal            = _mm256_and_pd(fscal,cutoff_mask);
 385
 386             /* Calculate temporary vectorial force */
 387             tx               = _mm256_mul_pd(fscal,dx30);
 388             ty               = _mm256_mul_pd(fscal,dy30);
 389             tz               = _mm256_mul_pd(fscal,dz30);
 390
 391             /* Update vectorial force */
 392             fix3             = _mm256_add_pd(fix3,tx);
 393             fiy3             = _mm256_add_pd(fiy3,ty);
 394             fiz3             = _mm256_add_pd(fiz3,tz);
 395
 396             fjx0             = _mm256_add_pd(fjx0,tx);
 397             fjy0             = _mm256_add_pd(fjy0,ty);
 398             fjz0             = _mm256_add_pd(fjz0,tz);
 399
 400             }
 401
 402             fjptrA             = f+j_coord_offsetA;
 403             fjptrB             = f+j_coord_offsetB;
 404             fjptrC             = f+j_coord_offsetC;
 405             fjptrD             = f+j_coord_offsetD;
 406
 407             gmx_mm256_decrement_1rvec_4ptr_swizzle_pd(fjptrA,fjptrB,fjptrC,fjptrD,fjx0,fjy0,fjz0);
 408
 409             /* Inner loop uses 141 flops */
 410         }
 411
 412         if(jidx<j_index_end)
 413         {
 414
 415             /* Get j neighbor index, and coordinate index */
 416             jnrlistA         = jjnr[jidx];
 417             jnrlistB         = jjnr[jidx+1];
 418             jnrlistC         = jjnr[jidx+2];
 419             jnrlistD         = jjnr[jidx+3];
 420             /* Sign of each element will be negative for non-real atoms.
 421              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
 422              * so use it as val = _mm_andnot_pd(mask,val) to clear dummy entries.
 423              */
 424             tmpmask0 = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
 425
 426             tmpmask1 = _mm_permute_ps(tmpmask0,_GMX_MM_PERMUTE(3,3,2,2));
 427             tmpmask0 = _mm_permute_ps(tmpmask0,_GMX_MM_PERMUTE(1,1,0,0));
 428             dummy_mask = _mm256_castps_pd(gmx_mm256_set_m128(tmpmask1,tmpmask0));
 429
 430             jnrA       = (jnrlistA>=0) ? jnrlistA : 0;
 431             jnrB       = (jnrlistB>=0) ? jnrlistB : 0;
 432             jnrC       = (jnrlistC>=0) ? jnrlistC : 0;
 433             jnrD       = (jnrlistD>=0) ? jnrlistD : 0;
 434             j_coord_offsetA  = DIM*jnrA;
 435             j_coord_offsetB  = DIM*jnrB;
 436             j_coord_offsetC  = DIM*jnrC;
 437             j_coord_offsetD  = DIM*jnrD;
 438
 439             /* load j atom coordinates */
 440             gmx_mm256_load_1rvec_4ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
 441                                                  x+j_coord_offsetC,x+j_coord_offsetD,
 442                                                  &jx0,&jy0,&jz0);
 443
 444             /* Calculate displacement vector */
 445             dx10             = _mm256_sub_pd(ix1,jx0);
 446             dy10             = _mm256_sub_pd(iy1,jy0);
 447             dz10             = _mm256_sub_pd(iz1,jz0);
 448             dx20             = _mm256_sub_pd(ix2,jx0);
 449             dy20             = _mm256_sub_pd(iy2,jy0);
 450             dz20             = _mm256_sub_pd(iz2,jz0);
 451             dx30             = _mm256_sub_pd(ix3,jx0);
 452             dy30             = _mm256_sub_pd(iy3,jy0);
 453             dz30             = _mm256_sub_pd(iz3,jz0);
 454
 455             /* Calculate squared distance and things based on it */
 456             rsq10            = gmx_mm256_calc_rsq_pd(dx10,dy10,dz10);
 457             rsq20            = gmx_mm256_calc_rsq_pd(dx20,dy20,dz20);
 458             rsq30            = gmx_mm256_calc_rsq_pd(dx30,dy30,dz30);
 459
 460             rinv10           = gmx_mm256_invsqrt_pd(rsq10);
 461             rinv20           = gmx_mm256_invsqrt_pd(rsq20);
 462             rinv30           = gmx_mm256_invsqrt_pd(rsq30);
 463
 464             rinvsq10         = _mm256_mul_pd(rinv10,rinv10);
 465             rinvsq20         = _mm256_mul_pd(rinv20,rinv20);
 466             rinvsq30         = _mm256_mul_pd(rinv30,rinv30);
 467
 468             /* Load parameters for j particles */
 469             jq0              = gmx_mm256_load_4real_swizzle_pd(charge+jnrA+0,charge+jnrB+0,
 470                                                                  charge+jnrC+0,charge+jnrD+0);
 471
 472             fjx0             = _mm256_setzero_pd();
 473             fjy0             = _mm256_setzero_pd();
 474             fjz0             = _mm256_setzero_pd();
 475
 476             /**************************
 477              * CALCULATE INTERACTIONS *
 478              **************************/
 479
 480             if (gmx_mm256_any_lt(rsq10,rcutoff2))
 481             {
 482
 483             r10              = _mm256_mul_pd(rsq10,rinv10);
 484             r10              = _mm256_andnot_pd(dummy_mask,r10);
 485
 486             /* Compute parameters for interactions between i and j atoms */
 487             qq10             = _mm256_mul_pd(iq1,jq0);
 488
 489             /* EWALD ELECTROSTATICS */
 490
 491             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 492             ewrt             = _mm256_mul_pd(r10,ewtabscale);
 493             ewitab           = _mm256_cvttpd_epi32(ewrt);
 494             eweps            = _mm256_sub_pd(ewrt,_mm256_round_pd(ewrt, _MM_FROUND_FLOOR));
 495             ewitab           = _mm_slli_epi32(ewitab,2);
 496             ewtabF           = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 497             ewtabD           = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 498             ewtabV           = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,2) );
 499             ewtabFn          = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,3) );
 500             GMX_MM256_FULLTRANSPOSE4_PD(ewtabF,ewtabD,ewtabV,ewtabFn);
 501             felec            = _mm256_add_pd(ewtabF,_mm256_mul_pd(eweps,ewtabD));
 502             velec            = _mm256_sub_pd(ewtabV,_mm256_mul_pd(_mm256_mul_pd(ewtabhalfspace,eweps),_mm256_add_pd(ewtabF,felec)));
 503             velec            = _mm256_mul_pd(qq10,_mm256_sub_pd(_mm256_sub_pd(rinv10,sh_ewald),velec));
 504             felec            = _mm256_mul_pd(_mm256_mul_pd(qq10,rinv10),_mm256_sub_pd(rinvsq10,felec));
 505
 506             cutoff_mask      = _mm256_cmp_pd(rsq10,rcutoff2,_CMP_LT_OQ);
 507
 508             /* Update potential sum for this i atom from the interaction with this j atom. */
 509             velec            = _mm256_and_pd(velec,cutoff_mask);
 510             velec            = _mm256_andnot_pd(dummy_mask,velec);
 511             velecsum         = _mm256_add_pd(velecsum,velec);
 512
 513             fscal            = felec;
 514
 515             fscal            = _mm256_and_pd(fscal,cutoff_mask);
 516
 517             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
 518
 519             /* Calculate temporary vectorial force */
 520             tx               = _mm256_mul_pd(fscal,dx10);
 521             ty               = _mm256_mul_pd(fscal,dy10);
 522             tz               = _mm256_mul_pd(fscal,dz10);
 523
 524             /* Update vectorial force */
 525             fix1             = _mm256_add_pd(fix1,tx);
 526             fiy1             = _mm256_add_pd(fiy1,ty);
 527             fiz1             = _mm256_add_pd(fiz1,tz);
 528
 529             fjx0             = _mm256_add_pd(fjx0,tx);
 530             fjy0             = _mm256_add_pd(fjy0,ty);
 531             fjz0             = _mm256_add_pd(fjz0,tz);
 532
 533             }
 534
 535             /**************************
 536              * CALCULATE INTERACTIONS *
 537              **************************/
 538
 539             if (gmx_mm256_any_lt(rsq20,rcutoff2))
 540             {
 541
 542             r20              = _mm256_mul_pd(rsq20,rinv20);
 543             r20              = _mm256_andnot_pd(dummy_mask,r20);
 544
 545             /* Compute parameters for interactions between i and j atoms */
 546             qq20             = _mm256_mul_pd(iq2,jq0);
 547
 548             /* EWALD ELECTROSTATICS */
 549
 550             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 551             ewrt             = _mm256_mul_pd(r20,ewtabscale);
 552             ewitab           = _mm256_cvttpd_epi32(ewrt);
 553             eweps            = _mm256_sub_pd(ewrt,_mm256_round_pd(ewrt, _MM_FROUND_FLOOR));
 554             ewitab           = _mm_slli_epi32(ewitab,2);
 555             ewtabF           = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 556             ewtabD           = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 557             ewtabV           = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,2) );
 558             ewtabFn          = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,3) );
 559             GMX_MM256_FULLTRANSPOSE4_PD(ewtabF,ewtabD,ewtabV,ewtabFn);
 560             felec            = _mm256_add_pd(ewtabF,_mm256_mul_pd(eweps,ewtabD));
 561             velec            = _mm256_sub_pd(ewtabV,_mm256_mul_pd(_mm256_mul_pd(ewtabhalfspace,eweps),_mm256_add_pd(ewtabF,felec)));
 562             velec            = _mm256_mul_pd(qq20,_mm256_sub_pd(_mm256_sub_pd(rinv20,sh_ewald),velec));
 563             felec            = _mm256_mul_pd(_mm256_mul_pd(qq20,rinv20),_mm256_sub_pd(rinvsq20,felec));
 564
 565             cutoff_mask      = _mm256_cmp_pd(rsq20,rcutoff2,_CMP_LT_OQ);
 566
 567             /* Update potential sum for this i atom from the interaction with this j atom. */
 568             velec            = _mm256_and_pd(velec,cutoff_mask);
 569             velec            = _mm256_andnot_pd(dummy_mask,velec);
 570             velecsum         = _mm256_add_pd(velecsum,velec);
 571
 572             fscal            = felec;
 573
 574             fscal            = _mm256_and_pd(fscal,cutoff_mask);
 575
 576             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
 577
 578             /* Calculate temporary vectorial force */
 579             tx               = _mm256_mul_pd(fscal,dx20);
 580             ty               = _mm256_mul_pd(fscal,dy20);
 581             tz               = _mm256_mul_pd(fscal,dz20);
 582
 583             /* Update vectorial force */
 584             fix2             = _mm256_add_pd(fix2,tx);
 585             fiy2             = _mm256_add_pd(fiy2,ty);
 586             fiz2             = _mm256_add_pd(fiz2,tz);
 587
 588             fjx0             = _mm256_add_pd(fjx0,tx);
 589             fjy0             = _mm256_add_pd(fjy0,ty);
 590             fjz0             = _mm256_add_pd(fjz0,tz);
 591
 592             }
 593
 594             /**************************
 595              * CALCULATE INTERACTIONS *
 596              **************************/
 597
 598             if (gmx_mm256_any_lt(rsq30,rcutoff2))
 599             {
 600
 601             r30              = _mm256_mul_pd(rsq30,rinv30);
 602             r30              = _mm256_andnot_pd(dummy_mask,r30);
 603
 604             /* Compute parameters for interactions between i and j atoms */
 605             qq30             = _mm256_mul_pd(iq3,jq0);
 606
 607             /* EWALD ELECTROSTATICS */
 608
 609             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 610             ewrt             = _mm256_mul_pd(r30,ewtabscale);
 611             ewitab           = _mm256_cvttpd_epi32(ewrt);
 612             eweps            = _mm256_sub_pd(ewrt,_mm256_round_pd(ewrt, _MM_FROUND_FLOOR));
 613             ewitab           = _mm_slli_epi32(ewitab,2);
 614             ewtabF           = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 615             ewtabD           = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 616             ewtabV           = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,2) );
 617             ewtabFn          = _mm256_load_pd( ewtab + _mm_extract_epi32(ewitab,3) );
 618             GMX_MM256_FULLTRANSPOSE4_PD(ewtabF,ewtabD,ewtabV,ewtabFn);
 619             felec            = _mm256_add_pd(ewtabF,_mm256_mul_pd(eweps,ewtabD));
 620             velec            = _mm256_sub_pd(ewtabV,_mm256_mul_pd(_mm256_mul_pd(ewtabhalfspace,eweps),_mm256_add_pd(ewtabF,felec)));
 621             velec            = _mm256_mul_pd(qq30,_mm256_sub_pd(_mm256_sub_pd(rinv30,sh_ewald),velec));
 622             felec            = _mm256_mul_pd(_mm256_mul_pd(qq30,rinv30),_mm256_sub_pd(rinvsq30,felec));
 623
 624             cutoff_mask      = _mm256_cmp_pd(rsq30,rcutoff2,_CMP_LT_OQ);
 625
 626             /* Update potential sum for this i atom from the interaction with this j atom. */
 627             velec            = _mm256_and_pd(velec,cutoff_mask);
 628             velec            = _mm256_andnot_pd(dummy_mask,velec);
 629             velecsum         = _mm256_add_pd(velecsum,velec);
 630
 631             fscal            = felec;
 632
 633             fscal            = _mm256_and_pd(fscal,cutoff_mask);
 634
 635             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
 636
 637             /* Calculate temporary vectorial force */
 638             tx               = _mm256_mul_pd(fscal,dx30);
 639             ty               = _mm256_mul_pd(fscal,dy30);
 640             tz               = _mm256_mul_pd(fscal,dz30);
 641
 642             /* Update vectorial force */
 643             fix3             = _mm256_add_pd(fix3,tx);
 644             fiy3             = _mm256_add_pd(fiy3,ty);
 645             fiz3             = _mm256_add_pd(fiz3,tz);
 646
 647             fjx0             = _mm256_add_pd(fjx0,tx);
 648             fjy0             = _mm256_add_pd(fjy0,ty);
 649             fjz0             = _mm256_add_pd(fjz0,tz);
 650
 651             }
 652
 653             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
 654             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
 655             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
 656             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
 657
 658             gmx_mm256_decrement_1rvec_4ptr_swizzle_pd(fjptrA,fjptrB,fjptrC,fjptrD,fjx0,fjy0,fjz0);
 659
 660             /* Inner loop uses 144 flops */
 661         }
 662
 663         /* End of innermost loop */
 664
 665         gmx_mm256_update_iforce_3atom_swizzle_pd(fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
 666                                                  f+i_coord_offset+DIM,fshift+i_shift_offset);
 667
 668         ggid                        = gid[iidx];
 669         /* Update potential energies */
 670         gmx_mm256_update_1pot_pd(velecsum,kernel_data->energygrp_elec+ggid);
 671
 672         /* Increment number of inner iterations */
 673         inneriter                  += j_index_end - j_index_start;
 674
 675         /* Outer loop uses 19 flops */
 676     }
 677
 678     /* Increment number of outer iterations */
 679     outeriter        += nri;
 680
 681     /* Update outer/inner flops */
 682
 683     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_W4_VF,outeriter*19 + inneriter*144);
 684 }
 685 /*
 686  * Gromacs nonbonded kernel:   nb_kernel_ElecEwSh_VdwNone_GeomW4P1_F_avx_256_double
 687  * Electrostatics interaction: Ewald
 688  * VdW interaction:            None
 689  * Geometry:                   Water4-Particle
 690  * Calculate force/pot:        Force
 691  */
 692 void
 693 nb_kernel_ElecEwSh_VdwNone_GeomW4P1_F_avx_256_double
 694                     (t_nblist                    * gmx_restrict       nlist,
 695                      rvec                        * gmx_restrict          xx,
 696                      rvec                        * gmx_restrict          ff,
 697                      t_forcerec                  * gmx_restrict          fr,
 698                      t_mdatoms                   * gmx_restrict     mdatoms,
 699                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
 700                      t_nrnb                      * gmx_restrict        nrnb)
 701 {
 702     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
 703      * just 0 for non-waters.
 704      * Suffixes A,B,C,D refer to j loop unrolling done with AVX, e.g. for the four different
 705      * jnr indices corresponding to data put in the four positions in the SIMD register.
 706      */
 707     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
 708     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
 709     int              jnrA,jnrB,jnrC,jnrD;
 710     int              jnrlistA,jnrlistB,jnrlistC,jnrlistD;
 711     int              jnrlistE,jnrlistF,jnrlistG,jnrlistH;
 712     int              j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
 713     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
 714     real             rcutoff_scalar;
 715     real             *shiftvec,*fshift,*x,*f;
 716     real             *fjptrA,*fjptrB,*fjptrC,*fjptrD;
 717     real             scratch[4*DIM];
 718     __m256d          tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
 719     real *           vdwioffsetptr1;
 720     __m256d          ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
 721     real *           vdwioffsetptr2;
 722     __m256d          ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
 723     real *           vdwioffsetptr3;
 724     __m256d          ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
 725     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
 726     __m256d          jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
 727     __m256d          dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
 728     __m256d          dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
 729     __m256d          dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
 730     __m256d          velec,felec,velecsum,facel,crf,krf,krf2;
 731     real             *charge;
 732     __m128i          ewitab;
 733     __m256d          ewtabscale,eweps,sh_ewald,ewrt,ewtabhalfspace,ewtabF,ewtabFn,ewtabD,ewtabV;
 734     __m256d          beta,beta2,beta3,zeta2,pmecorrF,pmecorrV,rinv3;
 735     real             *ewtab;
 736     __m256d          dummy_mask,cutoff_mask;
 737     __m128           tmpmask0,tmpmask1;
 738     __m256d          signbit = _mm256_castsi256_pd( _mm256_set1_epi32(0x80000000) );
 739     __m256d          one     = _mm256_set1_pd(1.0);
 740     __m256d          two     = _mm256_set1_pd(2.0);
 741     x                = xx[0];
 742     f                = ff[0];
 743
 744     nri              = nlist->nri;
 745     iinr             = nlist->iinr;
 746     jindex           = nlist->jindex;
 747     jjnr             = nlist->jjnr;
 748     shiftidx         = nlist->shift;
 749     gid              = nlist->gid;
 750     shiftvec         = fr->shift_vec[0];
 751     fshift           = fr->fshift[0];
 752     facel            = _mm256_set1_pd(fr->epsfac);
 753     charge           = mdatoms->chargeA;
 754
 755     sh_ewald         = _mm256_set1_pd(fr->ic->sh_ewald);
 756     beta             = _mm256_set1_pd(fr->ic->ewaldcoeff_q);
 757     beta2            = _mm256_mul_pd(beta,beta);
 758     beta3            = _mm256_mul_pd(beta,beta2);
 759
 760     ewtab            = fr->ic->tabq_coul_F;
 761     ewtabscale       = _mm256_set1_pd(fr->ic->tabq_scale);
 762     ewtabhalfspace   = _mm256_set1_pd(0.5/fr->ic->tabq_scale);
 763
 764     /* Setup water-specific parameters */
 765     inr              = nlist->iinr[0];
 766     iq1              = _mm256_mul_pd(facel,_mm256_set1_pd(charge[inr+1]));
 767     iq2              = _mm256_mul_pd(facel,_mm256_set1_pd(charge[inr+2]));
 768     iq3              = _mm256_mul_pd(facel,_mm256_set1_pd(charge[inr+3]));
 769
 770     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
 771     rcutoff_scalar   = fr->rcoulomb;
 772     rcutoff          = _mm256_set1_pd(rcutoff_scalar);
 773     rcutoff2         = _mm256_mul_pd(rcutoff,rcutoff);
 774
 775     /* Avoid stupid compiler warnings */
 776     jnrA = jnrB = jnrC = jnrD = 0;
 777     j_coord_offsetA = 0;
 778     j_coord_offsetB = 0;
 779     j_coord_offsetC = 0;
 780     j_coord_offsetD = 0;
 781
 782     outeriter        = 0;
 783     inneriter        = 0;
 784
 785     for(iidx=0;iidx<4*DIM;iidx++)
 786     {
 787         scratch[iidx] = 0.0;
 788     }
 789
 790     /* Start outer loop over neighborlists */
 791     for(iidx=0; iidx<nri; iidx++)
 792     {
 793         /* Load shift vector for this list */
 794         i_shift_offset   = DIM*shiftidx[iidx];
 795
 796         /* Load limits for loop over neighbors */
 797         j_index_start    = jindex[iidx];
 798         j_index_end      = jindex[iidx+1];
 799
 800         /* Get outer coordinate index */
 801         inr              = iinr[iidx];
 802         i_coord_offset   = DIM*inr;
 803
 804         /* Load i particle coords and add shift vector */
 805         gmx_mm256_load_shift_and_3rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset+DIM,
 806                                                     &ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
 807
 808         fix1             = _mm256_setzero_pd();
 809         fiy1             = _mm256_setzero_pd();
 810         fiz1             = _mm256_setzero_pd();
 811         fix2             = _mm256_setzero_pd();
 812         fiy2             = _mm256_setzero_pd();
 813         fiz2             = _mm256_setzero_pd();
 814         fix3             = _mm256_setzero_pd();
 815         fiy3             = _mm256_setzero_pd();
 816         fiz3             = _mm256_setzero_pd();
 817
 818         /* Start inner kernel loop */
 819         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
 820         {
 821
 822             /* Get j neighbor index, and coordinate index */
 823             jnrA             = jjnr[jidx];
 824             jnrB             = jjnr[jidx+1];
 825             jnrC             = jjnr[jidx+2];
 826             jnrD             = jjnr[jidx+3];
 827             j_coord_offsetA  = DIM*jnrA;
 828             j_coord_offsetB  = DIM*jnrB;
 829             j_coord_offsetC  = DIM*jnrC;
 830             j_coord_offsetD  = DIM*jnrD;
 831
 832             /* load j atom coordinates */
 833             gmx_mm256_load_1rvec_4ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
 834                                                  x+j_coord_offsetC,x+j_coord_offsetD,
 835                                                  &jx0,&jy0,&jz0);
 836
 837             /* Calculate displacement vector */
 838             dx10             = _mm256_sub_pd(ix1,jx0);
 839             dy10             = _mm256_sub_pd(iy1,jy0);
 840             dz10             = _mm256_sub_pd(iz1,jz0);
 841             dx20             = _mm256_sub_pd(ix2,jx0);
 842             dy20             = _mm256_sub_pd(iy2,jy0);
 843             dz20             = _mm256_sub_pd(iz2,jz0);
 844             dx30             = _mm256_sub_pd(ix3,jx0);
 845             dy30             = _mm256_sub_pd(iy3,jy0);
 846             dz30             = _mm256_sub_pd(iz3,jz0);
 847
 848             /* Calculate squared distance and things based on it */
 849             rsq10            = gmx_mm256_calc_rsq_pd(dx10,dy10,dz10);
 850             rsq20            = gmx_mm256_calc_rsq_pd(dx20,dy20,dz20);
 851             rsq30            = gmx_mm256_calc_rsq_pd(dx30,dy30,dz30);
 852
 853             rinv10           = gmx_mm256_invsqrt_pd(rsq10);
 854             rinv20           = gmx_mm256_invsqrt_pd(rsq20);
 855             rinv30           = gmx_mm256_invsqrt_pd(rsq30);
 856
 857             rinvsq10         = _mm256_mul_pd(rinv10,rinv10);
 858             rinvsq20         = _mm256_mul_pd(rinv20,rinv20);
 859             rinvsq30         = _mm256_mul_pd(rinv30,rinv30);
 860
 861             /* Load parameters for j particles */
 862             jq0              = gmx_mm256_load_4real_swizzle_pd(charge+jnrA+0,charge+jnrB+0,
 863                                                                  charge+jnrC+0,charge+jnrD+0);
 864
 865             fjx0             = _mm256_setzero_pd();
 866             fjy0             = _mm256_setzero_pd();
 867             fjz0             = _mm256_setzero_pd();
 868
 869             /**************************
 870              * CALCULATE INTERACTIONS *
 871              **************************/
 872
 873             if (gmx_mm256_any_lt(rsq10,rcutoff2))
 874             {
 875
 876             r10              = _mm256_mul_pd(rsq10,rinv10);
 877
 878             /* Compute parameters for interactions between i and j atoms */
 879             qq10             = _mm256_mul_pd(iq1,jq0);
 880
 881             /* EWALD ELECTROSTATICS */
 882
 883             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 884             ewrt             = _mm256_mul_pd(r10,ewtabscale);
 885             ewitab           = _mm256_cvttpd_epi32(ewrt);
 886             eweps            = _mm256_sub_pd(ewrt,_mm256_round_pd(ewrt, _MM_FROUND_FLOOR));
 887             gmx_mm256_load_4pair_swizzle_pd(ewtab + _mm_extract_epi32(ewitab,0),ewtab + _mm_extract_epi32(ewitab,1),
 888                                             ewtab + _mm_extract_epi32(ewitab,2),ewtab + _mm_extract_epi32(ewitab,3),
 889                                             &ewtabF,&ewtabFn);
 890             felec            = _mm256_add_pd(_mm256_mul_pd( _mm256_sub_pd(one,eweps),ewtabF),_mm256_mul_pd(eweps,ewtabFn));
 891             felec            = _mm256_mul_pd(_mm256_mul_pd(qq10,rinv10),_mm256_sub_pd(rinvsq10,felec));
 892
 893             cutoff_mask      = _mm256_cmp_pd(rsq10,rcutoff2,_CMP_LT_OQ);
 894
 895             fscal            = felec;
 896
 897             fscal            = _mm256_and_pd(fscal,cutoff_mask);
 898
 899             /* Calculate temporary vectorial force */
 900             tx               = _mm256_mul_pd(fscal,dx10);
 901             ty               = _mm256_mul_pd(fscal,dy10);
 902             tz               = _mm256_mul_pd(fscal,dz10);
 903
 904             /* Update vectorial force */
 905             fix1             = _mm256_add_pd(fix1,tx);
 906             fiy1             = _mm256_add_pd(fiy1,ty);
 907             fiz1             = _mm256_add_pd(fiz1,tz);
 908
 909             fjx0             = _mm256_add_pd(fjx0,tx);
 910             fjy0             = _mm256_add_pd(fjy0,ty);
 911             fjz0             = _mm256_add_pd(fjz0,tz);
 912
 913             }
 914
 915             /**************************
 916              * CALCULATE INTERACTIONS *
 917              **************************/
 918
 919             if (gmx_mm256_any_lt(rsq20,rcutoff2))
 920             {
 921
 922             r20              = _mm256_mul_pd(rsq20,rinv20);
 923
 924             /* Compute parameters for interactions between i and j atoms */
 925             qq20             = _mm256_mul_pd(iq2,jq0);
 926
 927             /* EWALD ELECTROSTATICS */
 928
 929             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 930             ewrt             = _mm256_mul_pd(r20,ewtabscale);
 931             ewitab           = _mm256_cvttpd_epi32(ewrt);
 932             eweps            = _mm256_sub_pd(ewrt,_mm256_round_pd(ewrt, _MM_FROUND_FLOOR));
 933             gmx_mm256_load_4pair_swizzle_pd(ewtab + _mm_extract_epi32(ewitab,0),ewtab + _mm_extract_epi32(ewitab,1),
 934                                             ewtab + _mm_extract_epi32(ewitab,2),ewtab + _mm_extract_epi32(ewitab,3),
 935                                             &ewtabF,&ewtabFn);
 936             felec            = _mm256_add_pd(_mm256_mul_pd( _mm256_sub_pd(one,eweps),ewtabF),_mm256_mul_pd(eweps,ewtabFn));
 937             felec            = _mm256_mul_pd(_mm256_mul_pd(qq20,rinv20),_mm256_sub_pd(rinvsq20,felec));
 938
 939             cutoff_mask      = _mm256_cmp_pd(rsq20,rcutoff2,_CMP_LT_OQ);
 940
 941             fscal            = felec;
 942
 943             fscal            = _mm256_and_pd(fscal,cutoff_mask);
 944
 945             /* Calculate temporary vectorial force */
 946             tx               = _mm256_mul_pd(fscal,dx20);
 947             ty               = _mm256_mul_pd(fscal,dy20);
 948             tz               = _mm256_mul_pd(fscal,dz20);
 949
 950             /* Update vectorial force */
 951             fix2             = _mm256_add_pd(fix2,tx);
 952             fiy2             = _mm256_add_pd(fiy2,ty);
 953             fiz2             = _mm256_add_pd(fiz2,tz);
 954
 955             fjx0             = _mm256_add_pd(fjx0,tx);
 956             fjy0             = _mm256_add_pd(fjy0,ty);
 957             fjz0             = _mm256_add_pd(fjz0,tz);
 958
 959             }
 960
 961             /**************************
 962              * CALCULATE INTERACTIONS *
 963              **************************/
 964
 965             if (gmx_mm256_any_lt(rsq30,rcutoff2))
 966             {
 967
 968             r30              = _mm256_mul_pd(rsq30,rinv30);
 969
 970             /* Compute parameters for interactions between i and j atoms */
 971             qq30             = _mm256_mul_pd(iq3,jq0);
 972
 973             /* EWALD ELECTROSTATICS */
 974
 975             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 976             ewrt             = _mm256_mul_pd(r30,ewtabscale);
 977             ewitab           = _mm256_cvttpd_epi32(ewrt);
 978             eweps            = _mm256_sub_pd(ewrt,_mm256_round_pd(ewrt, _MM_FROUND_FLOOR));
 979             gmx_mm256_load_4pair_swizzle_pd(ewtab + _mm_extract_epi32(ewitab,0),ewtab + _mm_extract_epi32(ewitab,1),
 980                                             ewtab + _mm_extract_epi32(ewitab,2),ewtab + _mm_extract_epi32(ewitab,3),
 981                                             &ewtabF,&ewtabFn);
 982             felec            = _mm256_add_pd(_mm256_mul_pd( _mm256_sub_pd(one,eweps),ewtabF),_mm256_mul_pd(eweps,ewtabFn));
 983             felec            = _mm256_mul_pd(_mm256_mul_pd(qq30,rinv30),_mm256_sub_pd(rinvsq30,felec));
 984
 985             cutoff_mask      = _mm256_cmp_pd(rsq30,rcutoff2,_CMP_LT_OQ);
 986
 987             fscal            = felec;
 988
 989             fscal            = _mm256_and_pd(fscal,cutoff_mask);
 990
 991             /* Calculate temporary vectorial force */
 992             tx               = _mm256_mul_pd(fscal,dx30);
 993             ty               = _mm256_mul_pd(fscal,dy30);
 994             tz               = _mm256_mul_pd(fscal,dz30);
 995
 996             /* Update vectorial force */
 997             fix3             = _mm256_add_pd(fix3,tx);
 998             fiy3             = _mm256_add_pd(fiy3,ty);
 999             fiz3             = _mm256_add_pd(fiz3,tz);
1000
1001             fjx0             = _mm256_add_pd(fjx0,tx);
1002             fjy0             = _mm256_add_pd(fjy0,ty);
1003             fjz0             = _mm256_add_pd(fjz0,tz);
1004
1005             }
1006
1007             fjptrA             = f+j_coord_offsetA;
1008             fjptrB             = f+j_coord_offsetB;
1009             fjptrC             = f+j_coord_offsetC;
1010             fjptrD             = f+j_coord_offsetD;
1011
1012             gmx_mm256_decrement_1rvec_4ptr_swizzle_pd(fjptrA,fjptrB,fjptrC,fjptrD,fjx0,fjy0,fjz0);
1013
1014             /* Inner loop uses 120 flops */
1015         }
1016
1017         if(jidx<j_index_end)
1018         {
1019
1020             /* Get j neighbor index, and coordinate index */
1021             jnrlistA         = jjnr[jidx];
1022             jnrlistB         = jjnr[jidx+1];
1023             jnrlistC         = jjnr[jidx+2];
1024             jnrlistD         = jjnr[jidx+3];
1025             /* Sign of each element will be negative for non-real atoms.
1026              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
1027              * so use it as val = _mm_andnot_pd(mask,val) to clear dummy entries.
1028              */
1029             tmpmask0 = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
1030
1031             tmpmask1 = _mm_permute_ps(tmpmask0,_GMX_MM_PERMUTE(3,3,2,2));
1032             tmpmask0 = _mm_permute_ps(tmpmask0,_GMX_MM_PERMUTE(1,1,0,0));
1033             dummy_mask = _mm256_castps_pd(gmx_mm256_set_m128(tmpmask1,tmpmask0));
1034
1035             jnrA       = (jnrlistA>=0) ? jnrlistA : 0;
1036             jnrB       = (jnrlistB>=0) ? jnrlistB : 0;
1037             jnrC       = (jnrlistC>=0) ? jnrlistC : 0;
1038             jnrD       = (jnrlistD>=0) ? jnrlistD : 0;
1039             j_coord_offsetA  = DIM*jnrA;
1040             j_coord_offsetB  = DIM*jnrB;
1041             j_coord_offsetC  = DIM*jnrC;
1042             j_coord_offsetD  = DIM*jnrD;
1043
1044             /* load j atom coordinates */
1045             gmx_mm256_load_1rvec_4ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
1046                                                  x+j_coord_offsetC,x+j_coord_offsetD,
1047                                                  &jx0,&jy0,&jz0);
1048
1049             /* Calculate displacement vector */
1050             dx10             = _mm256_sub_pd(ix1,jx0);
1051             dy10             = _mm256_sub_pd(iy1,jy0);
1052             dz10             = _mm256_sub_pd(iz1,jz0);
1053             dx20             = _mm256_sub_pd(ix2,jx0);
1054             dy20             = _mm256_sub_pd(iy2,jy0);
1055             dz20             = _mm256_sub_pd(iz2,jz0);
1056             dx30             = _mm256_sub_pd(ix3,jx0);
1057             dy30             = _mm256_sub_pd(iy3,jy0);
1058             dz30             = _mm256_sub_pd(iz3,jz0);
1059
1060             /* Calculate squared distance and things based on it */
1061             rsq10            = gmx_mm256_calc_rsq_pd(dx10,dy10,dz10);
1062             rsq20            = gmx_mm256_calc_rsq_pd(dx20,dy20,dz20);
1063             rsq30            = gmx_mm256_calc_rsq_pd(dx30,dy30,dz30);
1064
1065             rinv10           = gmx_mm256_invsqrt_pd(rsq10);
1066             rinv20           = gmx_mm256_invsqrt_pd(rsq20);
1067             rinv30           = gmx_mm256_invsqrt_pd(rsq30);
1068
1069             rinvsq10         = _mm256_mul_pd(rinv10,rinv10);
1070             rinvsq20         = _mm256_mul_pd(rinv20,rinv20);
1071             rinvsq30         = _mm256_mul_pd(rinv30,rinv30);
1072
1073             /* Load parameters for j particles */
1074             jq0              = gmx_mm256_load_4real_swizzle_pd(charge+jnrA+0,charge+jnrB+0,
1075                                                                  charge+jnrC+0,charge+jnrD+0);
1076
1077             fjx0             = _mm256_setzero_pd();
1078             fjy0             = _mm256_setzero_pd();
1079             fjz0             = _mm256_setzero_pd();
1080
1081             /**************************
1082              * CALCULATE INTERACTIONS *
1083              **************************/
1084
1085             if (gmx_mm256_any_lt(rsq10,rcutoff2))
1086             {
1087
1088             r10              = _mm256_mul_pd(rsq10,rinv10);
1089             r10              = _mm256_andnot_pd(dummy_mask,r10);
1090
1091             /* Compute parameters for interactions between i and j atoms */
1092             qq10             = _mm256_mul_pd(iq1,jq0);
1093
1094             /* EWALD ELECTROSTATICS */
1095
1096             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1097             ewrt             = _mm256_mul_pd(r10,ewtabscale);
1098             ewitab           = _mm256_cvttpd_epi32(ewrt);
1099             eweps            = _mm256_sub_pd(ewrt,_mm256_round_pd(ewrt, _MM_FROUND_FLOOR));
1100             gmx_mm256_load_4pair_swizzle_pd(ewtab + _mm_extract_epi32(ewitab,0),ewtab + _mm_extract_epi32(ewitab,1),
1101                                             ewtab + _mm_extract_epi32(ewitab,2),ewtab + _mm_extract_epi32(ewitab,3),
1102                                             &ewtabF,&ewtabFn);
1103             felec            = _mm256_add_pd(_mm256_mul_pd( _mm256_sub_pd(one,eweps),ewtabF),_mm256_mul_pd(eweps,ewtabFn));
1104             felec            = _mm256_mul_pd(_mm256_mul_pd(qq10,rinv10),_mm256_sub_pd(rinvsq10,felec));
1105
1106             cutoff_mask      = _mm256_cmp_pd(rsq10,rcutoff2,_CMP_LT_OQ);
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,dx10);
1116             ty               = _mm256_mul_pd(fscal,dy10);
1117             tz               = _mm256_mul_pd(fscal,dz10);
1118
1119             /* Update vectorial force */
1120             fix1             = _mm256_add_pd(fix1,tx);
1121             fiy1             = _mm256_add_pd(fiy1,ty);
1122             fiz1             = _mm256_add_pd(fiz1,tz);
1123
1124             fjx0             = _mm256_add_pd(fjx0,tx);
1125             fjy0             = _mm256_add_pd(fjy0,ty);
1126             fjz0             = _mm256_add_pd(fjz0,tz);
1127
1128             }
1129
1130             /**************************
1131              * CALCULATE INTERACTIONS *
1132              **************************/
1133
1134             if (gmx_mm256_any_lt(rsq20,rcutoff2))
1135             {
1136
1137             r20              = _mm256_mul_pd(rsq20,rinv20);
1138             r20              = _mm256_andnot_pd(dummy_mask,r20);
1139
1140             /* Compute parameters for interactions between i and j atoms */
1141             qq20             = _mm256_mul_pd(iq2,jq0);
1142
1143             /* EWALD ELECTROSTATICS */
1144
1145             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1146             ewrt             = _mm256_mul_pd(r20,ewtabscale);
1147             ewitab           = _mm256_cvttpd_epi32(ewrt);
1148             eweps            = _mm256_sub_pd(ewrt,_mm256_round_pd(ewrt, _MM_FROUND_FLOOR));
1149             gmx_mm256_load_4pair_swizzle_pd(ewtab + _mm_extract_epi32(ewitab,0),ewtab + _mm_extract_epi32(ewitab,1),
1150                                             ewtab + _mm_extract_epi32(ewitab,2),ewtab + _mm_extract_epi32(ewitab,3),
1151                                             &ewtabF,&ewtabFn);
1152             felec            = _mm256_add_pd(_mm256_mul_pd( _mm256_sub_pd(one,eweps),ewtabF),_mm256_mul_pd(eweps,ewtabFn));
1153             felec            = _mm256_mul_pd(_mm256_mul_pd(qq20,rinv20),_mm256_sub_pd(rinvsq20,felec));
1154
1155             cutoff_mask      = _mm256_cmp_pd(rsq20,rcutoff2,_CMP_LT_OQ);
1156
1157             fscal            = felec;
1158
1159             fscal            = _mm256_and_pd(fscal,cutoff_mask);
1160
1161             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
1162
1163             /* Calculate temporary vectorial force */
1164             tx               = _mm256_mul_pd(fscal,dx20);
1165             ty               = _mm256_mul_pd(fscal,dy20);
1166             tz               = _mm256_mul_pd(fscal,dz20);
1167
1168             /* Update vectorial force */
1169             fix2             = _mm256_add_pd(fix2,tx);
1170             fiy2             = _mm256_add_pd(fiy2,ty);
1171             fiz2             = _mm256_add_pd(fiz2,tz);
1172
1173             fjx0             = _mm256_add_pd(fjx0,tx);
1174             fjy0             = _mm256_add_pd(fjy0,ty);
1175             fjz0             = _mm256_add_pd(fjz0,tz);
1176
1177             }
1178
1179             /**************************
1180              * CALCULATE INTERACTIONS *
1181              **************************/
1182
1183             if (gmx_mm256_any_lt(rsq30,rcutoff2))
1184             {
1185
1186             r30              = _mm256_mul_pd(rsq30,rinv30);
1187             r30              = _mm256_andnot_pd(dummy_mask,r30);
1188
1189             /* Compute parameters for interactions between i and j atoms */
1190             qq30             = _mm256_mul_pd(iq3,jq0);
1191
1192             /* EWALD ELECTROSTATICS */
1193
1194             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1195             ewrt             = _mm256_mul_pd(r30,ewtabscale);
1196             ewitab           = _mm256_cvttpd_epi32(ewrt);
1197             eweps            = _mm256_sub_pd(ewrt,_mm256_round_pd(ewrt, _MM_FROUND_FLOOR));
1198             gmx_mm256_load_4pair_swizzle_pd(ewtab + _mm_extract_epi32(ewitab,0),ewtab + _mm_extract_epi32(ewitab,1),
1199                                             ewtab + _mm_extract_epi32(ewitab,2),ewtab + _mm_extract_epi32(ewitab,3),
1200                                             &ewtabF,&ewtabFn);
1201             felec            = _mm256_add_pd(_mm256_mul_pd( _mm256_sub_pd(one,eweps),ewtabF),_mm256_mul_pd(eweps,ewtabFn));
1202             felec            = _mm256_mul_pd(_mm256_mul_pd(qq30,rinv30),_mm256_sub_pd(rinvsq30,felec));
1203
1204             cutoff_mask      = _mm256_cmp_pd(rsq30,rcutoff2,_CMP_LT_OQ);
1205
1206             fscal            = felec;
1207
1208             fscal            = _mm256_and_pd(fscal,cutoff_mask);
1209
1210             fscal            = _mm256_andnot_pd(dummy_mask,fscal);
1211
1212             /* Calculate temporary vectorial force */
1213             tx               = _mm256_mul_pd(fscal,dx30);
1214             ty               = _mm256_mul_pd(fscal,dy30);
1215             tz               = _mm256_mul_pd(fscal,dz30);
1216
1217             /* Update vectorial force */
1218             fix3             = _mm256_add_pd(fix3,tx);
1219             fiy3             = _mm256_add_pd(fiy3,ty);
1220             fiz3             = _mm256_add_pd(fiz3,tz);
1221
1222             fjx0             = _mm256_add_pd(fjx0,tx);
1223             fjy0             = _mm256_add_pd(fjy0,ty);
1224             fjz0             = _mm256_add_pd(fjz0,tz);
1225
1226             }
1227
1228             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
1229             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
1230             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
1231             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
1232
1233             gmx_mm256_decrement_1rvec_4ptr_swizzle_pd(fjptrA,fjptrB,fjptrC,fjptrD,fjx0,fjy0,fjz0);
1234
1235             /* Inner loop uses 123 flops */
1236         }
1237
1238         /* End of innermost loop */
1239
1240         gmx_mm256_update_iforce_3atom_swizzle_pd(fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
1241                                                  f+i_coord_offset+DIM,fshift+i_shift_offset);
1242
1243         /* Increment number of inner iterations */
1244         inneriter                  += j_index_end - j_index_start;
1245
1246         /* Outer loop uses 18 flops */
1247     }
1248
1249     /* Increment number of outer iterations */
1250     outeriter        += nri;
1251
1252     /* Update outer/inner flops */
1253
1254     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_W4_F,outeriter*18 + inneriter*123);
1255 }