src/gromacs/gmxlib/nonbonded/nb_kernel_avx_128_fma_double/nb_kernel_ElecEwSw_VdwLJSw_GeomW4P1_avx_128_fma_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_128_fma_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_128_fma_double.h"
  48 #include "kernelutil_x86_avx_128_fma_double.h"
  49
  50 /*
  51  * Gromacs nonbonded kernel:   nb_kernel_ElecEwSw_VdwLJSw_GeomW4P1_VF_avx_128_fma_double
  52  * Electrostatics interaction: Ewald
  53  * VdW interaction:            LennardJones
  54  * Geometry:                   Water4-Particle
  55  * Calculate force/pot:        PotentialAndForce
  56  */
  57 void
  58 nb_kernel_ElecEwSw_VdwLJSw_GeomW4P1_VF_avx_128_fma_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 refer to j loop unrolling done with SSE double precision, e.g. for the two 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;
  75     int              j_coord_offsetA,j_coord_offsetB;
  76     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
  77     real             rcutoff_scalar;
  78     real             *shiftvec,*fshift,*x,*f;
  79     __m128d          tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
  80     int              vdwioffset0;
  81     __m128d          ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
  82     int              vdwioffset1;
  83     __m128d          ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
  84     int              vdwioffset2;
  85     __m128d          ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
  86     int              vdwioffset3;
  87     __m128d          ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
  88     int              vdwjidx0A,vdwjidx0B;
  89     __m128d          jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
  90     __m128d          dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
  91     __m128d          dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
  92     __m128d          dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
  93     __m128d          dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
  94     __m128d          velec,felec,velecsum,facel,crf,krf,krf2;
  95     real             *charge;
  96     int              nvdwtype;
  97     __m128d          rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
  98     int              *vdwtype;
  99     real             *vdwparam;
 100     __m128d          one_sixth   = _mm_set1_pd(1.0/6.0);
 101     __m128d          one_twelfth = _mm_set1_pd(1.0/12.0);
 102     __m128i          ewitab;
 103     __m128d          ewtabscale,eweps,twoeweps,sh_ewald,ewrt,ewtabhalfspace,ewtabF,ewtabFn,ewtabD,ewtabV;
 104     real             *ewtab;
 105     __m128d          rswitch,swV3,swV4,swV5,swF2,swF3,swF4,d,d2,sw,dsw;
 106     real             rswitch_scalar,d_scalar;
 107     __m128d          dummy_mask,cutoff_mask;
 108     __m128d          signbit   = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
 109     __m128d          one     = _mm_set1_pd(1.0);
 110     __m128d          two     = _mm_set1_pd(2.0);
 111     x                = xx[0];
 112     f                = ff[0];
 113
 114     nri              = nlist->nri;
 115     iinr             = nlist->iinr;
 116     jindex           = nlist->jindex;
 117     jjnr             = nlist->jjnr;
 118     shiftidx         = nlist->shift;
 119     gid              = nlist->gid;
 120     shiftvec         = fr->shift_vec[0];
 121     fshift           = fr->fshift[0];
 122     facel            = _mm_set1_pd(fr->epsfac);
 123     charge           = mdatoms->chargeA;
 124     nvdwtype         = fr->ntype;
 125     vdwparam         = fr->nbfp;
 126     vdwtype          = mdatoms->typeA;
 127
 128     sh_ewald         = _mm_set1_pd(fr->ic->sh_ewald);
 129     ewtab            = fr->ic->tabq_coul_FDV0;
 130     ewtabscale       = _mm_set1_pd(fr->ic->tabq_scale);
 131     ewtabhalfspace   = _mm_set1_pd(0.5/fr->ic->tabq_scale);
 132
 133     /* Setup water-specific parameters */
 134     inr              = nlist->iinr[0];
 135     iq1              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+1]));
 136     iq2              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+2]));
 137     iq3              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+3]));
 138     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
 139
 140     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
 141     rcutoff_scalar   = fr->rcoulomb;
 142     rcutoff          = _mm_set1_pd(rcutoff_scalar);
 143     rcutoff2         = _mm_mul_pd(rcutoff,rcutoff);
 144
 145     rswitch_scalar   = fr->rcoulomb_switch;
 146     rswitch          = _mm_set1_pd(rswitch_scalar);
 147     /* Setup switch parameters */
 148     d_scalar         = rcutoff_scalar-rswitch_scalar;
 149     d                = _mm_set1_pd(d_scalar);
 150     swV3             = _mm_set1_pd(-10.0/(d_scalar*d_scalar*d_scalar));
 151     swV4             = _mm_set1_pd( 15.0/(d_scalar*d_scalar*d_scalar*d_scalar));
 152     swV5             = _mm_set1_pd( -6.0/(d_scalar*d_scalar*d_scalar*d_scalar*d_scalar));
 153     swF2             = _mm_set1_pd(-30.0/(d_scalar*d_scalar*d_scalar));
 154     swF3             = _mm_set1_pd( 60.0/(d_scalar*d_scalar*d_scalar*d_scalar));
 155     swF4             = _mm_set1_pd(-30.0/(d_scalar*d_scalar*d_scalar*d_scalar*d_scalar));
 156
 157     /* Avoid stupid compiler warnings */
 158     jnrA = jnrB = 0;
 159     j_coord_offsetA = 0;
 160     j_coord_offsetB = 0;
 161
 162     outeriter        = 0;
 163     inneriter        = 0;
 164
 165     /* Start outer loop over neighborlists */
 166     for(iidx=0; iidx<nri; iidx++)
 167     {
 168         /* Load shift vector for this list */
 169         i_shift_offset   = DIM*shiftidx[iidx];
 170
 171         /* Load limits for loop over neighbors */
 172         j_index_start    = jindex[iidx];
 173         j_index_end      = jindex[iidx+1];
 174
 175         /* Get outer coordinate index */
 176         inr              = iinr[iidx];
 177         i_coord_offset   = DIM*inr;
 178
 179         /* Load i particle coords and add shift vector */
 180         gmx_mm_load_shift_and_4rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,
 181                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
 182
 183         fix0             = _mm_setzero_pd();
 184         fiy0             = _mm_setzero_pd();
 185         fiz0             = _mm_setzero_pd();
 186         fix1             = _mm_setzero_pd();
 187         fiy1             = _mm_setzero_pd();
 188         fiz1             = _mm_setzero_pd();
 189         fix2             = _mm_setzero_pd();
 190         fiy2             = _mm_setzero_pd();
 191         fiz2             = _mm_setzero_pd();
 192         fix3             = _mm_setzero_pd();
 193         fiy3             = _mm_setzero_pd();
 194         fiz3             = _mm_setzero_pd();
 195
 196         /* Reset potential sums */
 197         velecsum         = _mm_setzero_pd();
 198         vvdwsum          = _mm_setzero_pd();
 199
 200         /* Start inner kernel loop */
 201         for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
 202         {
 203
 204             /* Get j neighbor index, and coordinate index */
 205             jnrA             = jjnr[jidx];
 206             jnrB             = jjnr[jidx+1];
 207             j_coord_offsetA  = DIM*jnrA;
 208             j_coord_offsetB  = DIM*jnrB;
 209
 210             /* load j atom coordinates */
 211             gmx_mm_load_1rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
 212                                               &jx0,&jy0,&jz0);
 213
 214             /* Calculate displacement vector */
 215             dx00             = _mm_sub_pd(ix0,jx0);
 216             dy00             = _mm_sub_pd(iy0,jy0);
 217             dz00             = _mm_sub_pd(iz0,jz0);
 218             dx10             = _mm_sub_pd(ix1,jx0);
 219             dy10             = _mm_sub_pd(iy1,jy0);
 220             dz10             = _mm_sub_pd(iz1,jz0);
 221             dx20             = _mm_sub_pd(ix2,jx0);
 222             dy20             = _mm_sub_pd(iy2,jy0);
 223             dz20             = _mm_sub_pd(iz2,jz0);
 224             dx30             = _mm_sub_pd(ix3,jx0);
 225             dy30             = _mm_sub_pd(iy3,jy0);
 226             dz30             = _mm_sub_pd(iz3,jz0);
 227
 228             /* Calculate squared distance and things based on it */
 229             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
 230             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
 231             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
 232             rsq30            = gmx_mm_calc_rsq_pd(dx30,dy30,dz30);
 233
 234             rinv00           = gmx_mm_invsqrt_pd(rsq00);
 235             rinv10           = gmx_mm_invsqrt_pd(rsq10);
 236             rinv20           = gmx_mm_invsqrt_pd(rsq20);
 237             rinv30           = gmx_mm_invsqrt_pd(rsq30);
 238
 239             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
 240             rinvsq10         = _mm_mul_pd(rinv10,rinv10);
 241             rinvsq20         = _mm_mul_pd(rinv20,rinv20);
 242             rinvsq30         = _mm_mul_pd(rinv30,rinv30);
 243
 244             /* Load parameters for j particles */
 245             jq0              = gmx_mm_load_2real_swizzle_pd(charge+jnrA+0,charge+jnrB+0);
 246             vdwjidx0A        = 2*vdwtype[jnrA+0];
 247             vdwjidx0B        = 2*vdwtype[jnrB+0];
 248
 249             fjx0             = _mm_setzero_pd();
 250             fjy0             = _mm_setzero_pd();
 251             fjz0             = _mm_setzero_pd();
 252
 253             /**************************
 254              * CALCULATE INTERACTIONS *
 255              **************************/
 256
 257             if (gmx_mm_any_lt(rsq00,rcutoff2))
 258             {
 259
 260             r00              = _mm_mul_pd(rsq00,rinv00);
 261
 262             /* Compute parameters for interactions between i and j atoms */
 263             gmx_mm_load_2pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,
 264                                          vdwparam+vdwioffset0+vdwjidx0B,&c6_00,&c12_00);
 265
 266             /* LENNARD-JONES DISPERSION/REPULSION */
 267
 268             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
 269             vvdw6            = _mm_mul_pd(c6_00,rinvsix);
 270             vvdw12           = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
 271             vvdw             = _mm_msub_pd( vvdw12,one_twelfth, _mm_mul_pd(vvdw6,one_sixth) );
 272             fvdw             = _mm_mul_pd(_mm_sub_pd(vvdw12,vvdw6),rinvsq00);
 273
 274             d                = _mm_sub_pd(r00,rswitch);
 275             d                = _mm_max_pd(d,_mm_setzero_pd());
 276             d2               = _mm_mul_pd(d,d);
 277             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
 278
 279             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
 280
 281             /* Evaluate switch function */
 282             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 283             fvdw             = _mm_msub_pd( fvdw,sw , _mm_mul_pd(rinv00,_mm_mul_pd(vvdw,dsw)) );
 284             vvdw             = _mm_mul_pd(vvdw,sw);
 285             cutoff_mask      = _mm_cmplt_pd(rsq00,rcutoff2);
 286
 287             /* Update potential sum for this i atom from the interaction with this j atom. */
 288             vvdw             = _mm_and_pd(vvdw,cutoff_mask);
 289             vvdwsum          = _mm_add_pd(vvdwsum,vvdw);
 290
 291             fscal            = fvdw;
 292
 293             fscal            = _mm_and_pd(fscal,cutoff_mask);
 294
 295             /* Update vectorial force */
 296             fix0             = _mm_macc_pd(dx00,fscal,fix0);
 297             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
 298             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
 299
 300             fjx0             = _mm_macc_pd(dx00,fscal,fjx0);
 301             fjy0             = _mm_macc_pd(dy00,fscal,fjy0);
 302             fjz0             = _mm_macc_pd(dz00,fscal,fjz0);
 303
 304             }
 305
 306             /**************************
 307              * CALCULATE INTERACTIONS *
 308              **************************/
 309
 310             if (gmx_mm_any_lt(rsq10,rcutoff2))
 311             {
 312
 313             r10              = _mm_mul_pd(rsq10,rinv10);
 314
 315             /* Compute parameters for interactions between i and j atoms */
 316             qq10             = _mm_mul_pd(iq1,jq0);
 317
 318             /* EWALD ELECTROSTATICS */
 319
 320             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 321             ewrt             = _mm_mul_pd(r10,ewtabscale);
 322             ewitab           = _mm_cvttpd_epi32(ewrt);
 323 #ifdef __XOP__
 324             eweps            = _mm_frcz_pd(ewrt);
 325 #else
 326             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 327 #endif
 328             twoeweps         = _mm_add_pd(eweps,eweps);
 329             ewitab           = _mm_slli_epi32(ewitab,2);
 330             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 331             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 332             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 333             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 334             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
 335             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 336             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 337             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 338             velec            = _mm_mul_pd(qq10,_mm_sub_pd(rinv10,velec));
 339             felec            = _mm_mul_pd(_mm_mul_pd(qq10,rinv10),_mm_sub_pd(rinvsq10,felec));
 340
 341             d                = _mm_sub_pd(r10,rswitch);
 342             d                = _mm_max_pd(d,_mm_setzero_pd());
 343             d2               = _mm_mul_pd(d,d);
 344             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
 345
 346             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
 347
 348             /* Evaluate switch function */
 349             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 350             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv10,_mm_mul_pd(velec,dsw)) );
 351             velec            = _mm_mul_pd(velec,sw);
 352             cutoff_mask      = _mm_cmplt_pd(rsq10,rcutoff2);
 353
 354             /* Update potential sum for this i atom from the interaction with this j atom. */
 355             velec            = _mm_and_pd(velec,cutoff_mask);
 356             velecsum         = _mm_add_pd(velecsum,velec);
 357
 358             fscal            = felec;
 359
 360             fscal            = _mm_and_pd(fscal,cutoff_mask);
 361
 362             /* Update vectorial force */
 363             fix1             = _mm_macc_pd(dx10,fscal,fix1);
 364             fiy1             = _mm_macc_pd(dy10,fscal,fiy1);
 365             fiz1             = _mm_macc_pd(dz10,fscal,fiz1);
 366
 367             fjx0             = _mm_macc_pd(dx10,fscal,fjx0);
 368             fjy0             = _mm_macc_pd(dy10,fscal,fjy0);
 369             fjz0             = _mm_macc_pd(dz10,fscal,fjz0);
 370
 371             }
 372
 373             /**************************
 374              * CALCULATE INTERACTIONS *
 375              **************************/
 376
 377             if (gmx_mm_any_lt(rsq20,rcutoff2))
 378             {
 379
 380             r20              = _mm_mul_pd(rsq20,rinv20);
 381
 382             /* Compute parameters for interactions between i and j atoms */
 383             qq20             = _mm_mul_pd(iq2,jq0);
 384
 385             /* EWALD ELECTROSTATICS */
 386
 387             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 388             ewrt             = _mm_mul_pd(r20,ewtabscale);
 389             ewitab           = _mm_cvttpd_epi32(ewrt);
 390 #ifdef __XOP__
 391             eweps            = _mm_frcz_pd(ewrt);
 392 #else
 393             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 394 #endif
 395             twoeweps         = _mm_add_pd(eweps,eweps);
 396             ewitab           = _mm_slli_epi32(ewitab,2);
 397             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 398             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 399             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 400             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 401             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
 402             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 403             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 404             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 405             velec            = _mm_mul_pd(qq20,_mm_sub_pd(rinv20,velec));
 406             felec            = _mm_mul_pd(_mm_mul_pd(qq20,rinv20),_mm_sub_pd(rinvsq20,felec));
 407
 408             d                = _mm_sub_pd(r20,rswitch);
 409             d                = _mm_max_pd(d,_mm_setzero_pd());
 410             d2               = _mm_mul_pd(d,d);
 411             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
 412
 413             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
 414
 415             /* Evaluate switch function */
 416             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 417             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv20,_mm_mul_pd(velec,dsw)) );
 418             velec            = _mm_mul_pd(velec,sw);
 419             cutoff_mask      = _mm_cmplt_pd(rsq20,rcutoff2);
 420
 421             /* Update potential sum for this i atom from the interaction with this j atom. */
 422             velec            = _mm_and_pd(velec,cutoff_mask);
 423             velecsum         = _mm_add_pd(velecsum,velec);
 424
 425             fscal            = felec;
 426
 427             fscal            = _mm_and_pd(fscal,cutoff_mask);
 428
 429             /* Update vectorial force */
 430             fix2             = _mm_macc_pd(dx20,fscal,fix2);
 431             fiy2             = _mm_macc_pd(dy20,fscal,fiy2);
 432             fiz2             = _mm_macc_pd(dz20,fscal,fiz2);
 433
 434             fjx0             = _mm_macc_pd(dx20,fscal,fjx0);
 435             fjy0             = _mm_macc_pd(dy20,fscal,fjy0);
 436             fjz0             = _mm_macc_pd(dz20,fscal,fjz0);
 437
 438             }
 439
 440             /**************************
 441              * CALCULATE INTERACTIONS *
 442              **************************/
 443
 444             if (gmx_mm_any_lt(rsq30,rcutoff2))
 445             {
 446
 447             r30              = _mm_mul_pd(rsq30,rinv30);
 448
 449             /* Compute parameters for interactions between i and j atoms */
 450             qq30             = _mm_mul_pd(iq3,jq0);
 451
 452             /* EWALD ELECTROSTATICS */
 453
 454             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 455             ewrt             = _mm_mul_pd(r30,ewtabscale);
 456             ewitab           = _mm_cvttpd_epi32(ewrt);
 457 #ifdef __XOP__
 458             eweps            = _mm_frcz_pd(ewrt);
 459 #else
 460             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 461 #endif
 462             twoeweps         = _mm_add_pd(eweps,eweps);
 463             ewitab           = _mm_slli_epi32(ewitab,2);
 464             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 465             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 466             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 467             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 468             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
 469             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 470             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 471             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 472             velec            = _mm_mul_pd(qq30,_mm_sub_pd(rinv30,velec));
 473             felec            = _mm_mul_pd(_mm_mul_pd(qq30,rinv30),_mm_sub_pd(rinvsq30,felec));
 474
 475             d                = _mm_sub_pd(r30,rswitch);
 476             d                = _mm_max_pd(d,_mm_setzero_pd());
 477             d2               = _mm_mul_pd(d,d);
 478             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
 479
 480             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
 481
 482             /* Evaluate switch function */
 483             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 484             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv30,_mm_mul_pd(velec,dsw)) );
 485             velec            = _mm_mul_pd(velec,sw);
 486             cutoff_mask      = _mm_cmplt_pd(rsq30,rcutoff2);
 487
 488             /* Update potential sum for this i atom from the interaction with this j atom. */
 489             velec            = _mm_and_pd(velec,cutoff_mask);
 490             velecsum         = _mm_add_pd(velecsum,velec);
 491
 492             fscal            = felec;
 493
 494             fscal            = _mm_and_pd(fscal,cutoff_mask);
 495
 496             /* Update vectorial force */
 497             fix3             = _mm_macc_pd(dx30,fscal,fix3);
 498             fiy3             = _mm_macc_pd(dy30,fscal,fiy3);
 499             fiz3             = _mm_macc_pd(dz30,fscal,fiz3);
 500
 501             fjx0             = _mm_macc_pd(dx30,fscal,fjx0);
 502             fjy0             = _mm_macc_pd(dy30,fscal,fjy0);
 503             fjz0             = _mm_macc_pd(dz30,fscal,fjz0);
 504
 505             }
 506
 507             gmx_mm_decrement_1rvec_2ptr_swizzle_pd(f+j_coord_offsetA,f+j_coord_offsetB,fjx0,fjy0,fjz0);
 508
 509             /* Inner loop uses 269 flops */
 510         }
 511
 512         if(jidx<j_index_end)
 513         {
 514
 515             jnrA             = jjnr[jidx];
 516             j_coord_offsetA  = DIM*jnrA;
 517
 518             /* load j atom coordinates */
 519             gmx_mm_load_1rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
 520                                               &jx0,&jy0,&jz0);
 521
 522             /* Calculate displacement vector */
 523             dx00             = _mm_sub_pd(ix0,jx0);
 524             dy00             = _mm_sub_pd(iy0,jy0);
 525             dz00             = _mm_sub_pd(iz0,jz0);
 526             dx10             = _mm_sub_pd(ix1,jx0);
 527             dy10             = _mm_sub_pd(iy1,jy0);
 528             dz10             = _mm_sub_pd(iz1,jz0);
 529             dx20             = _mm_sub_pd(ix2,jx0);
 530             dy20             = _mm_sub_pd(iy2,jy0);
 531             dz20             = _mm_sub_pd(iz2,jz0);
 532             dx30             = _mm_sub_pd(ix3,jx0);
 533             dy30             = _mm_sub_pd(iy3,jy0);
 534             dz30             = _mm_sub_pd(iz3,jz0);
 535
 536             /* Calculate squared distance and things based on it */
 537             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
 538             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
 539             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
 540             rsq30            = gmx_mm_calc_rsq_pd(dx30,dy30,dz30);
 541
 542             rinv00           = gmx_mm_invsqrt_pd(rsq00);
 543             rinv10           = gmx_mm_invsqrt_pd(rsq10);
 544             rinv20           = gmx_mm_invsqrt_pd(rsq20);
 545             rinv30           = gmx_mm_invsqrt_pd(rsq30);
 546
 547             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
 548             rinvsq10         = _mm_mul_pd(rinv10,rinv10);
 549             rinvsq20         = _mm_mul_pd(rinv20,rinv20);
 550             rinvsq30         = _mm_mul_pd(rinv30,rinv30);
 551
 552             /* Load parameters for j particles */
 553             jq0              = _mm_load_sd(charge+jnrA+0);
 554             vdwjidx0A        = 2*vdwtype[jnrA+0];
 555
 556             fjx0             = _mm_setzero_pd();
 557             fjy0             = _mm_setzero_pd();
 558             fjz0             = _mm_setzero_pd();
 559
 560             /**************************
 561              * CALCULATE INTERACTIONS *
 562              **************************/
 563
 564             if (gmx_mm_any_lt(rsq00,rcutoff2))
 565             {
 566
 567             r00              = _mm_mul_pd(rsq00,rinv00);
 568
 569             /* Compute parameters for interactions between i and j atoms */
 570             gmx_mm_load_1pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,&c6_00,&c12_00);
 571
 572             /* LENNARD-JONES DISPERSION/REPULSION */
 573
 574             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
 575             vvdw6            = _mm_mul_pd(c6_00,rinvsix);
 576             vvdw12           = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
 577             vvdw             = _mm_msub_pd( vvdw12,one_twelfth, _mm_mul_pd(vvdw6,one_sixth) );
 578             fvdw             = _mm_mul_pd(_mm_sub_pd(vvdw12,vvdw6),rinvsq00);
 579
 580             d                = _mm_sub_pd(r00,rswitch);
 581             d                = _mm_max_pd(d,_mm_setzero_pd());
 582             d2               = _mm_mul_pd(d,d);
 583             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
 584
 585             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
 586
 587             /* Evaluate switch function */
 588             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 589             fvdw             = _mm_msub_pd( fvdw,sw , _mm_mul_pd(rinv00,_mm_mul_pd(vvdw,dsw)) );
 590             vvdw             = _mm_mul_pd(vvdw,sw);
 591             cutoff_mask      = _mm_cmplt_pd(rsq00,rcutoff2);
 592
 593             /* Update potential sum for this i atom from the interaction with this j atom. */
 594             vvdw             = _mm_and_pd(vvdw,cutoff_mask);
 595             vvdw             = _mm_unpacklo_pd(vvdw,_mm_setzero_pd());
 596             vvdwsum          = _mm_add_pd(vvdwsum,vvdw);
 597
 598             fscal            = fvdw;
 599
 600             fscal            = _mm_and_pd(fscal,cutoff_mask);
 601
 602             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
 603
 604             /* Update vectorial force */
 605             fix0             = _mm_macc_pd(dx00,fscal,fix0);
 606             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
 607             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
 608
 609             fjx0             = _mm_macc_pd(dx00,fscal,fjx0);
 610             fjy0             = _mm_macc_pd(dy00,fscal,fjy0);
 611             fjz0             = _mm_macc_pd(dz00,fscal,fjz0);
 612
 613             }
 614
 615             /**************************
 616              * CALCULATE INTERACTIONS *
 617              **************************/
 618
 619             if (gmx_mm_any_lt(rsq10,rcutoff2))
 620             {
 621
 622             r10              = _mm_mul_pd(rsq10,rinv10);
 623
 624             /* Compute parameters for interactions between i and j atoms */
 625             qq10             = _mm_mul_pd(iq1,jq0);
 626
 627             /* EWALD ELECTROSTATICS */
 628
 629             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 630             ewrt             = _mm_mul_pd(r10,ewtabscale);
 631             ewitab           = _mm_cvttpd_epi32(ewrt);
 632 #ifdef __XOP__
 633             eweps            = _mm_frcz_pd(ewrt);
 634 #else
 635             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 636 #endif
 637             twoeweps         = _mm_add_pd(eweps,eweps);
 638             ewitab           = _mm_slli_epi32(ewitab,2);
 639             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 640             ewtabD           = _mm_setzero_pd();
 641             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 642             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 643             ewtabFn          = _mm_setzero_pd();
 644             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 645             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 646             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 647             velec            = _mm_mul_pd(qq10,_mm_sub_pd(rinv10,velec));
 648             felec            = _mm_mul_pd(_mm_mul_pd(qq10,rinv10),_mm_sub_pd(rinvsq10,felec));
 649
 650             d                = _mm_sub_pd(r10,rswitch);
 651             d                = _mm_max_pd(d,_mm_setzero_pd());
 652             d2               = _mm_mul_pd(d,d);
 653             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
 654
 655             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
 656
 657             /* Evaluate switch function */
 658             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 659             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv10,_mm_mul_pd(velec,dsw)) );
 660             velec            = _mm_mul_pd(velec,sw);
 661             cutoff_mask      = _mm_cmplt_pd(rsq10,rcutoff2);
 662
 663             /* Update potential sum for this i atom from the interaction with this j atom. */
 664             velec            = _mm_and_pd(velec,cutoff_mask);
 665             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
 666             velecsum         = _mm_add_pd(velecsum,velec);
 667
 668             fscal            = felec;
 669
 670             fscal            = _mm_and_pd(fscal,cutoff_mask);
 671
 672             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
 673
 674             /* Update vectorial force */
 675             fix1             = _mm_macc_pd(dx10,fscal,fix1);
 676             fiy1             = _mm_macc_pd(dy10,fscal,fiy1);
 677             fiz1             = _mm_macc_pd(dz10,fscal,fiz1);
 678
 679             fjx0             = _mm_macc_pd(dx10,fscal,fjx0);
 680             fjy0             = _mm_macc_pd(dy10,fscal,fjy0);
 681             fjz0             = _mm_macc_pd(dz10,fscal,fjz0);
 682
 683             }
 684
 685             /**************************
 686              * CALCULATE INTERACTIONS *
 687              **************************/
 688
 689             if (gmx_mm_any_lt(rsq20,rcutoff2))
 690             {
 691
 692             r20              = _mm_mul_pd(rsq20,rinv20);
 693
 694             /* Compute parameters for interactions between i and j atoms */
 695             qq20             = _mm_mul_pd(iq2,jq0);
 696
 697             /* EWALD ELECTROSTATICS */
 698
 699             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 700             ewrt             = _mm_mul_pd(r20,ewtabscale);
 701             ewitab           = _mm_cvttpd_epi32(ewrt);
 702 #ifdef __XOP__
 703             eweps            = _mm_frcz_pd(ewrt);
 704 #else
 705             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 706 #endif
 707             twoeweps         = _mm_add_pd(eweps,eweps);
 708             ewitab           = _mm_slli_epi32(ewitab,2);
 709             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 710             ewtabD           = _mm_setzero_pd();
 711             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 712             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 713             ewtabFn          = _mm_setzero_pd();
 714             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 715             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 716             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 717             velec            = _mm_mul_pd(qq20,_mm_sub_pd(rinv20,velec));
 718             felec            = _mm_mul_pd(_mm_mul_pd(qq20,rinv20),_mm_sub_pd(rinvsq20,felec));
 719
 720             d                = _mm_sub_pd(r20,rswitch);
 721             d                = _mm_max_pd(d,_mm_setzero_pd());
 722             d2               = _mm_mul_pd(d,d);
 723             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
 724
 725             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
 726
 727             /* Evaluate switch function */
 728             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 729             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv20,_mm_mul_pd(velec,dsw)) );
 730             velec            = _mm_mul_pd(velec,sw);
 731             cutoff_mask      = _mm_cmplt_pd(rsq20,rcutoff2);
 732
 733             /* Update potential sum for this i atom from the interaction with this j atom. */
 734             velec            = _mm_and_pd(velec,cutoff_mask);
 735             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
 736             velecsum         = _mm_add_pd(velecsum,velec);
 737
 738             fscal            = felec;
 739
 740             fscal            = _mm_and_pd(fscal,cutoff_mask);
 741
 742             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
 743
 744             /* Update vectorial force */
 745             fix2             = _mm_macc_pd(dx20,fscal,fix2);
 746             fiy2             = _mm_macc_pd(dy20,fscal,fiy2);
 747             fiz2             = _mm_macc_pd(dz20,fscal,fiz2);
 748
 749             fjx0             = _mm_macc_pd(dx20,fscal,fjx0);
 750             fjy0             = _mm_macc_pd(dy20,fscal,fjy0);
 751             fjz0             = _mm_macc_pd(dz20,fscal,fjz0);
 752
 753             }
 754
 755             /**************************
 756              * CALCULATE INTERACTIONS *
 757              **************************/
 758
 759             if (gmx_mm_any_lt(rsq30,rcutoff2))
 760             {
 761
 762             r30              = _mm_mul_pd(rsq30,rinv30);
 763
 764             /* Compute parameters for interactions between i and j atoms */
 765             qq30             = _mm_mul_pd(iq3,jq0);
 766
 767             /* EWALD ELECTROSTATICS */
 768
 769             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 770             ewrt             = _mm_mul_pd(r30,ewtabscale);
 771             ewitab           = _mm_cvttpd_epi32(ewrt);
 772 #ifdef __XOP__
 773             eweps            = _mm_frcz_pd(ewrt);
 774 #else
 775             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 776 #endif
 777             twoeweps         = _mm_add_pd(eweps,eweps);
 778             ewitab           = _mm_slli_epi32(ewitab,2);
 779             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 780             ewtabD           = _mm_setzero_pd();
 781             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 782             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 783             ewtabFn          = _mm_setzero_pd();
 784             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 785             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 786             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 787             velec            = _mm_mul_pd(qq30,_mm_sub_pd(rinv30,velec));
 788             felec            = _mm_mul_pd(_mm_mul_pd(qq30,rinv30),_mm_sub_pd(rinvsq30,felec));
 789
 790             d                = _mm_sub_pd(r30,rswitch);
 791             d                = _mm_max_pd(d,_mm_setzero_pd());
 792             d2               = _mm_mul_pd(d,d);
 793             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
 794
 795             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
 796
 797             /* Evaluate switch function */
 798             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 799             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv30,_mm_mul_pd(velec,dsw)) );
 800             velec            = _mm_mul_pd(velec,sw);
 801             cutoff_mask      = _mm_cmplt_pd(rsq30,rcutoff2);
 802
 803             /* Update potential sum for this i atom from the interaction with this j atom. */
 804             velec            = _mm_and_pd(velec,cutoff_mask);
 805             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
 806             velecsum         = _mm_add_pd(velecsum,velec);
 807
 808             fscal            = felec;
 809
 810             fscal            = _mm_and_pd(fscal,cutoff_mask);
 811
 812             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
 813
 814             /* Update vectorial force */
 815             fix3             = _mm_macc_pd(dx30,fscal,fix3);
 816             fiy3             = _mm_macc_pd(dy30,fscal,fiy3);
 817             fiz3             = _mm_macc_pd(dz30,fscal,fiz3);
 818
 819             fjx0             = _mm_macc_pd(dx30,fscal,fjx0);
 820             fjy0             = _mm_macc_pd(dy30,fscal,fjy0);
 821             fjz0             = _mm_macc_pd(dz30,fscal,fjz0);
 822
 823             }
 824
 825             gmx_mm_decrement_1rvec_1ptr_swizzle_pd(f+j_coord_offsetA,fjx0,fjy0,fjz0);
 826
 827             /* Inner loop uses 269 flops */
 828         }
 829
 830         /* End of innermost loop */
 831
 832         gmx_mm_update_iforce_4atom_swizzle_pd(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
 833                                               f+i_coord_offset,fshift+i_shift_offset);
 834
 835         ggid                        = gid[iidx];
 836         /* Update potential energies */
 837         gmx_mm_update_1pot_pd(velecsum,kernel_data->energygrp_elec+ggid);
 838         gmx_mm_update_1pot_pd(vvdwsum,kernel_data->energygrp_vdw+ggid);
 839
 840         /* Increment number of inner iterations */
 841         inneriter                  += j_index_end - j_index_start;
 842
 843         /* Outer loop uses 26 flops */
 844     }
 845
 846     /* Increment number of outer iterations */
 847     outeriter        += nri;
 848
 849     /* Update outer/inner flops */
 850
 851     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4_VF,outeriter*26 + inneriter*269);
 852 }
 853 /*
 854  * Gromacs nonbonded kernel:   nb_kernel_ElecEwSw_VdwLJSw_GeomW4P1_F_avx_128_fma_double
 855  * Electrostatics interaction: Ewald
 856  * VdW interaction:            LennardJones
 857  * Geometry:                   Water4-Particle
 858  * Calculate force/pot:        Force
 859  */
 860 void
 861 nb_kernel_ElecEwSw_VdwLJSw_GeomW4P1_F_avx_128_fma_double
 862                     (t_nblist                    * gmx_restrict       nlist,
 863                      rvec                        * gmx_restrict          xx,
 864                      rvec                        * gmx_restrict          ff,
 865                      t_forcerec                  * gmx_restrict          fr,
 866                      t_mdatoms                   * gmx_restrict     mdatoms,
 867                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
 868                      t_nrnb                      * gmx_restrict        nrnb)
 869 {
 870     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
 871      * just 0 for non-waters.
 872      * Suffixes A,B refer to j loop unrolling done with SSE double precision, e.g. for the two different
 873      * jnr indices corresponding to data put in the four positions in the SIMD register.
 874      */
 875     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
 876     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
 877     int              jnrA,jnrB;
 878     int              j_coord_offsetA,j_coord_offsetB;
 879     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
 880     real             rcutoff_scalar;
 881     real             *shiftvec,*fshift,*x,*f;
 882     __m128d          tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
 883     int              vdwioffset0;
 884     __m128d          ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
 885     int              vdwioffset1;
 886     __m128d          ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
 887     int              vdwioffset2;
 888     __m128d          ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
 889     int              vdwioffset3;
 890     __m128d          ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
 891     int              vdwjidx0A,vdwjidx0B;
 892     __m128d          jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
 893     __m128d          dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
 894     __m128d          dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
 895     __m128d          dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
 896     __m128d          dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
 897     __m128d          velec,felec,velecsum,facel,crf,krf,krf2;
 898     real             *charge;
 899     int              nvdwtype;
 900     __m128d          rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
 901     int              *vdwtype;
 902     real             *vdwparam;
 903     __m128d          one_sixth   = _mm_set1_pd(1.0/6.0);
 904     __m128d          one_twelfth = _mm_set1_pd(1.0/12.0);
 905     __m128i          ewitab;
 906     __m128d          ewtabscale,eweps,twoeweps,sh_ewald,ewrt,ewtabhalfspace,ewtabF,ewtabFn,ewtabD,ewtabV;
 907     real             *ewtab;
 908     __m128d          rswitch,swV3,swV4,swV5,swF2,swF3,swF4,d,d2,sw,dsw;
 909     real             rswitch_scalar,d_scalar;
 910     __m128d          dummy_mask,cutoff_mask;
 911     __m128d          signbit   = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
 912     __m128d          one     = _mm_set1_pd(1.0);
 913     __m128d          two     = _mm_set1_pd(2.0);
 914     x                = xx[0];
 915     f                = ff[0];
 916
 917     nri              = nlist->nri;
 918     iinr             = nlist->iinr;
 919     jindex           = nlist->jindex;
 920     jjnr             = nlist->jjnr;
 921     shiftidx         = nlist->shift;
 922     gid              = nlist->gid;
 923     shiftvec         = fr->shift_vec[0];
 924     fshift           = fr->fshift[0];
 925     facel            = _mm_set1_pd(fr->epsfac);
 926     charge           = mdatoms->chargeA;
 927     nvdwtype         = fr->ntype;
 928     vdwparam         = fr->nbfp;
 929     vdwtype          = mdatoms->typeA;
 930
 931     sh_ewald         = _mm_set1_pd(fr->ic->sh_ewald);
 932     ewtab            = fr->ic->tabq_coul_FDV0;
 933     ewtabscale       = _mm_set1_pd(fr->ic->tabq_scale);
 934     ewtabhalfspace   = _mm_set1_pd(0.5/fr->ic->tabq_scale);
 935
 936     /* Setup water-specific parameters */
 937     inr              = nlist->iinr[0];
 938     iq1              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+1]));
 939     iq2              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+2]));
 940     iq3              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+3]));
 941     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
 942
 943     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
 944     rcutoff_scalar   = fr->rcoulomb;
 945     rcutoff          = _mm_set1_pd(rcutoff_scalar);
 946     rcutoff2         = _mm_mul_pd(rcutoff,rcutoff);
 947
 948     rswitch_scalar   = fr->rcoulomb_switch;
 949     rswitch          = _mm_set1_pd(rswitch_scalar);
 950     /* Setup switch parameters */
 951     d_scalar         = rcutoff_scalar-rswitch_scalar;
 952     d                = _mm_set1_pd(d_scalar);
 953     swV3             = _mm_set1_pd(-10.0/(d_scalar*d_scalar*d_scalar));
 954     swV4             = _mm_set1_pd( 15.0/(d_scalar*d_scalar*d_scalar*d_scalar));
 955     swV5             = _mm_set1_pd( -6.0/(d_scalar*d_scalar*d_scalar*d_scalar*d_scalar));
 956     swF2             = _mm_set1_pd(-30.0/(d_scalar*d_scalar*d_scalar));
 957     swF3             = _mm_set1_pd( 60.0/(d_scalar*d_scalar*d_scalar*d_scalar));
 958     swF4             = _mm_set1_pd(-30.0/(d_scalar*d_scalar*d_scalar*d_scalar*d_scalar));
 959
 960     /* Avoid stupid compiler warnings */
 961     jnrA = jnrB = 0;
 962     j_coord_offsetA = 0;
 963     j_coord_offsetB = 0;
 964
 965     outeriter        = 0;
 966     inneriter        = 0;
 967
 968     /* Start outer loop over neighborlists */
 969     for(iidx=0; iidx<nri; iidx++)
 970     {
 971         /* Load shift vector for this list */
 972         i_shift_offset   = DIM*shiftidx[iidx];
 973
 974         /* Load limits for loop over neighbors */
 975         j_index_start    = jindex[iidx];
 976         j_index_end      = jindex[iidx+1];
 977
 978         /* Get outer coordinate index */
 979         inr              = iinr[iidx];
 980         i_coord_offset   = DIM*inr;
 981
 982         /* Load i particle coords and add shift vector */
 983         gmx_mm_load_shift_and_4rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,
 984                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
 985
 986         fix0             = _mm_setzero_pd();
 987         fiy0             = _mm_setzero_pd();
 988         fiz0             = _mm_setzero_pd();
 989         fix1             = _mm_setzero_pd();
 990         fiy1             = _mm_setzero_pd();
 991         fiz1             = _mm_setzero_pd();
 992         fix2             = _mm_setzero_pd();
 993         fiy2             = _mm_setzero_pd();
 994         fiz2             = _mm_setzero_pd();
 995         fix3             = _mm_setzero_pd();
 996         fiy3             = _mm_setzero_pd();
 997         fiz3             = _mm_setzero_pd();
 998
 999         /* Start inner kernel loop */
1000         for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
1001         {
1002
1003             /* Get j neighbor index, and coordinate index */
1004             jnrA             = jjnr[jidx];
1005             jnrB             = jjnr[jidx+1];
1006             j_coord_offsetA  = DIM*jnrA;
1007             j_coord_offsetB  = DIM*jnrB;
1008
1009             /* load j atom coordinates */
1010             gmx_mm_load_1rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
1011                                               &jx0,&jy0,&jz0);
1012
1013             /* Calculate displacement vector */
1014             dx00             = _mm_sub_pd(ix0,jx0);
1015             dy00             = _mm_sub_pd(iy0,jy0);
1016             dz00             = _mm_sub_pd(iz0,jz0);
1017             dx10             = _mm_sub_pd(ix1,jx0);
1018             dy10             = _mm_sub_pd(iy1,jy0);
1019             dz10             = _mm_sub_pd(iz1,jz0);
1020             dx20             = _mm_sub_pd(ix2,jx0);
1021             dy20             = _mm_sub_pd(iy2,jy0);
1022             dz20             = _mm_sub_pd(iz2,jz0);
1023             dx30             = _mm_sub_pd(ix3,jx0);
1024             dy30             = _mm_sub_pd(iy3,jy0);
1025             dz30             = _mm_sub_pd(iz3,jz0);
1026
1027             /* Calculate squared distance and things based on it */
1028             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
1029             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
1030             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
1031             rsq30            = gmx_mm_calc_rsq_pd(dx30,dy30,dz30);
1032
1033             rinv00           = gmx_mm_invsqrt_pd(rsq00);
1034             rinv10           = gmx_mm_invsqrt_pd(rsq10);
1035             rinv20           = gmx_mm_invsqrt_pd(rsq20);
1036             rinv30           = gmx_mm_invsqrt_pd(rsq30);
1037
1038             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
1039             rinvsq10         = _mm_mul_pd(rinv10,rinv10);
1040             rinvsq20         = _mm_mul_pd(rinv20,rinv20);
1041             rinvsq30         = _mm_mul_pd(rinv30,rinv30);
1042
1043             /* Load parameters for j particles */
1044             jq0              = gmx_mm_load_2real_swizzle_pd(charge+jnrA+0,charge+jnrB+0);
1045             vdwjidx0A        = 2*vdwtype[jnrA+0];
1046             vdwjidx0B        = 2*vdwtype[jnrB+0];
1047
1048             fjx0             = _mm_setzero_pd();
1049             fjy0             = _mm_setzero_pd();
1050             fjz0             = _mm_setzero_pd();
1051
1052             /**************************
1053              * CALCULATE INTERACTIONS *
1054              **************************/
1055
1056             if (gmx_mm_any_lt(rsq00,rcutoff2))
1057             {
1058
1059             r00              = _mm_mul_pd(rsq00,rinv00);
1060
1061             /* Compute parameters for interactions between i and j atoms */
1062             gmx_mm_load_2pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,
1063                                          vdwparam+vdwioffset0+vdwjidx0B,&c6_00,&c12_00);
1064
1065             /* LENNARD-JONES DISPERSION/REPULSION */
1066
1067             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
1068             vvdw6            = _mm_mul_pd(c6_00,rinvsix);
1069             vvdw12           = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
1070             vvdw             = _mm_msub_pd( vvdw12,one_twelfth, _mm_mul_pd(vvdw6,one_sixth) );
1071             fvdw             = _mm_mul_pd(_mm_sub_pd(vvdw12,vvdw6),rinvsq00);
1072
1073             d                = _mm_sub_pd(r00,rswitch);
1074             d                = _mm_max_pd(d,_mm_setzero_pd());
1075             d2               = _mm_mul_pd(d,d);
1076             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
1077
1078             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
1079
1080             /* Evaluate switch function */
1081             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1082             fvdw             = _mm_msub_pd( fvdw,sw , _mm_mul_pd(rinv00,_mm_mul_pd(vvdw,dsw)) );
1083             cutoff_mask      = _mm_cmplt_pd(rsq00,rcutoff2);
1084
1085             fscal            = fvdw;
1086
1087             fscal            = _mm_and_pd(fscal,cutoff_mask);
1088
1089             /* Update vectorial force */
1090             fix0             = _mm_macc_pd(dx00,fscal,fix0);
1091             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
1092             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
1093
1094             fjx0             = _mm_macc_pd(dx00,fscal,fjx0);
1095             fjy0             = _mm_macc_pd(dy00,fscal,fjy0);
1096             fjz0             = _mm_macc_pd(dz00,fscal,fjz0);
1097
1098             }
1099
1100             /**************************
1101              * CALCULATE INTERACTIONS *
1102              **************************/
1103
1104             if (gmx_mm_any_lt(rsq10,rcutoff2))
1105             {
1106
1107             r10              = _mm_mul_pd(rsq10,rinv10);
1108
1109             /* Compute parameters for interactions between i and j atoms */
1110             qq10             = _mm_mul_pd(iq1,jq0);
1111
1112             /* EWALD ELECTROSTATICS */
1113
1114             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1115             ewrt             = _mm_mul_pd(r10,ewtabscale);
1116             ewitab           = _mm_cvttpd_epi32(ewrt);
1117 #ifdef __XOP__
1118             eweps            = _mm_frcz_pd(ewrt);
1119 #else
1120             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1121 #endif
1122             twoeweps         = _mm_add_pd(eweps,eweps);
1123             ewitab           = _mm_slli_epi32(ewitab,2);
1124             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
1125             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
1126             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1127             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
1128             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
1129             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1130             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
1131             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
1132             velec            = _mm_mul_pd(qq10,_mm_sub_pd(rinv10,velec));
1133             felec            = _mm_mul_pd(_mm_mul_pd(qq10,rinv10),_mm_sub_pd(rinvsq10,felec));
1134
1135             d                = _mm_sub_pd(r10,rswitch);
1136             d                = _mm_max_pd(d,_mm_setzero_pd());
1137             d2               = _mm_mul_pd(d,d);
1138             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
1139
1140             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
1141
1142             /* Evaluate switch function */
1143             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1144             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv10,_mm_mul_pd(velec,dsw)) );
1145             cutoff_mask      = _mm_cmplt_pd(rsq10,rcutoff2);
1146
1147             fscal            = felec;
1148
1149             fscal            = _mm_and_pd(fscal,cutoff_mask);
1150
1151             /* Update vectorial force */
1152             fix1             = _mm_macc_pd(dx10,fscal,fix1);
1153             fiy1             = _mm_macc_pd(dy10,fscal,fiy1);
1154             fiz1             = _mm_macc_pd(dz10,fscal,fiz1);
1155
1156             fjx0             = _mm_macc_pd(dx10,fscal,fjx0);
1157             fjy0             = _mm_macc_pd(dy10,fscal,fjy0);
1158             fjz0             = _mm_macc_pd(dz10,fscal,fjz0);
1159
1160             }
1161
1162             /**************************
1163              * CALCULATE INTERACTIONS *
1164              **************************/
1165
1166             if (gmx_mm_any_lt(rsq20,rcutoff2))
1167             {
1168
1169             r20              = _mm_mul_pd(rsq20,rinv20);
1170
1171             /* Compute parameters for interactions between i and j atoms */
1172             qq20             = _mm_mul_pd(iq2,jq0);
1173
1174             /* EWALD ELECTROSTATICS */
1175
1176             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1177             ewrt             = _mm_mul_pd(r20,ewtabscale);
1178             ewitab           = _mm_cvttpd_epi32(ewrt);
1179 #ifdef __XOP__
1180             eweps            = _mm_frcz_pd(ewrt);
1181 #else
1182             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1183 #endif
1184             twoeweps         = _mm_add_pd(eweps,eweps);
1185             ewitab           = _mm_slli_epi32(ewitab,2);
1186             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
1187             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
1188             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1189             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
1190             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
1191             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1192             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
1193             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
1194             velec            = _mm_mul_pd(qq20,_mm_sub_pd(rinv20,velec));
1195             felec            = _mm_mul_pd(_mm_mul_pd(qq20,rinv20),_mm_sub_pd(rinvsq20,felec));
1196
1197             d                = _mm_sub_pd(r20,rswitch);
1198             d                = _mm_max_pd(d,_mm_setzero_pd());
1199             d2               = _mm_mul_pd(d,d);
1200             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
1201
1202             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
1203
1204             /* Evaluate switch function */
1205             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1206             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv20,_mm_mul_pd(velec,dsw)) );
1207             cutoff_mask      = _mm_cmplt_pd(rsq20,rcutoff2);
1208
1209             fscal            = felec;
1210
1211             fscal            = _mm_and_pd(fscal,cutoff_mask);
1212
1213             /* Update vectorial force */
1214             fix2             = _mm_macc_pd(dx20,fscal,fix2);
1215             fiy2             = _mm_macc_pd(dy20,fscal,fiy2);
1216             fiz2             = _mm_macc_pd(dz20,fscal,fiz2);
1217
1218             fjx0             = _mm_macc_pd(dx20,fscal,fjx0);
1219             fjy0             = _mm_macc_pd(dy20,fscal,fjy0);
1220             fjz0             = _mm_macc_pd(dz20,fscal,fjz0);
1221
1222             }
1223
1224             /**************************
1225              * CALCULATE INTERACTIONS *
1226              **************************/
1227
1228             if (gmx_mm_any_lt(rsq30,rcutoff2))
1229             {
1230
1231             r30              = _mm_mul_pd(rsq30,rinv30);
1232
1233             /* Compute parameters for interactions between i and j atoms */
1234             qq30             = _mm_mul_pd(iq3,jq0);
1235
1236             /* EWALD ELECTROSTATICS */
1237
1238             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1239             ewrt             = _mm_mul_pd(r30,ewtabscale);
1240             ewitab           = _mm_cvttpd_epi32(ewrt);
1241 #ifdef __XOP__
1242             eweps            = _mm_frcz_pd(ewrt);
1243 #else
1244             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1245 #endif
1246             twoeweps         = _mm_add_pd(eweps,eweps);
1247             ewitab           = _mm_slli_epi32(ewitab,2);
1248             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
1249             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
1250             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1251             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
1252             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
1253             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1254             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
1255             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
1256             velec            = _mm_mul_pd(qq30,_mm_sub_pd(rinv30,velec));
1257             felec            = _mm_mul_pd(_mm_mul_pd(qq30,rinv30),_mm_sub_pd(rinvsq30,felec));
1258
1259             d                = _mm_sub_pd(r30,rswitch);
1260             d                = _mm_max_pd(d,_mm_setzero_pd());
1261             d2               = _mm_mul_pd(d,d);
1262             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
1263
1264             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
1265
1266             /* Evaluate switch function */
1267             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1268             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv30,_mm_mul_pd(velec,dsw)) );
1269             cutoff_mask      = _mm_cmplt_pd(rsq30,rcutoff2);
1270
1271             fscal            = felec;
1272
1273             fscal            = _mm_and_pd(fscal,cutoff_mask);
1274
1275             /* Update vectorial force */
1276             fix3             = _mm_macc_pd(dx30,fscal,fix3);
1277             fiy3             = _mm_macc_pd(dy30,fscal,fiy3);
1278             fiz3             = _mm_macc_pd(dz30,fscal,fiz3);
1279
1280             fjx0             = _mm_macc_pd(dx30,fscal,fjx0);
1281             fjy0             = _mm_macc_pd(dy30,fscal,fjy0);
1282             fjz0             = _mm_macc_pd(dz30,fscal,fjz0);
1283
1284             }
1285
1286             gmx_mm_decrement_1rvec_2ptr_swizzle_pd(f+j_coord_offsetA,f+j_coord_offsetB,fjx0,fjy0,fjz0);
1287
1288             /* Inner loop uses 257 flops */
1289         }
1290
1291         if(jidx<j_index_end)
1292         {
1293
1294             jnrA             = jjnr[jidx];
1295             j_coord_offsetA  = DIM*jnrA;
1296
1297             /* load j atom coordinates */
1298             gmx_mm_load_1rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
1299                                               &jx0,&jy0,&jz0);
1300
1301             /* Calculate displacement vector */
1302             dx00             = _mm_sub_pd(ix0,jx0);
1303             dy00             = _mm_sub_pd(iy0,jy0);
1304             dz00             = _mm_sub_pd(iz0,jz0);
1305             dx10             = _mm_sub_pd(ix1,jx0);
1306             dy10             = _mm_sub_pd(iy1,jy0);
1307             dz10             = _mm_sub_pd(iz1,jz0);
1308             dx20             = _mm_sub_pd(ix2,jx0);
1309             dy20             = _mm_sub_pd(iy2,jy0);
1310             dz20             = _mm_sub_pd(iz2,jz0);
1311             dx30             = _mm_sub_pd(ix3,jx0);
1312             dy30             = _mm_sub_pd(iy3,jy0);
1313             dz30             = _mm_sub_pd(iz3,jz0);
1314
1315             /* Calculate squared distance and things based on it */
1316             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
1317             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
1318             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
1319             rsq30            = gmx_mm_calc_rsq_pd(dx30,dy30,dz30);
1320
1321             rinv00           = gmx_mm_invsqrt_pd(rsq00);
1322             rinv10           = gmx_mm_invsqrt_pd(rsq10);
1323             rinv20           = gmx_mm_invsqrt_pd(rsq20);
1324             rinv30           = gmx_mm_invsqrt_pd(rsq30);
1325
1326             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
1327             rinvsq10         = _mm_mul_pd(rinv10,rinv10);
1328             rinvsq20         = _mm_mul_pd(rinv20,rinv20);
1329             rinvsq30         = _mm_mul_pd(rinv30,rinv30);
1330
1331             /* Load parameters for j particles */
1332             jq0              = _mm_load_sd(charge+jnrA+0);
1333             vdwjidx0A        = 2*vdwtype[jnrA+0];
1334
1335             fjx0             = _mm_setzero_pd();
1336             fjy0             = _mm_setzero_pd();
1337             fjz0             = _mm_setzero_pd();
1338
1339             /**************************
1340              * CALCULATE INTERACTIONS *
1341              **************************/
1342
1343             if (gmx_mm_any_lt(rsq00,rcutoff2))
1344             {
1345
1346             r00              = _mm_mul_pd(rsq00,rinv00);
1347
1348             /* Compute parameters for interactions between i and j atoms */
1349             gmx_mm_load_1pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,&c6_00,&c12_00);
1350
1351             /* LENNARD-JONES DISPERSION/REPULSION */
1352
1353             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
1354             vvdw6            = _mm_mul_pd(c6_00,rinvsix);
1355             vvdw12           = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
1356             vvdw             = _mm_msub_pd( vvdw12,one_twelfth, _mm_mul_pd(vvdw6,one_sixth) );
1357             fvdw             = _mm_mul_pd(_mm_sub_pd(vvdw12,vvdw6),rinvsq00);
1358
1359             d                = _mm_sub_pd(r00,rswitch);
1360             d                = _mm_max_pd(d,_mm_setzero_pd());
1361             d2               = _mm_mul_pd(d,d);
1362             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
1363
1364             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
1365
1366             /* Evaluate switch function */
1367             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1368             fvdw             = _mm_msub_pd( fvdw,sw , _mm_mul_pd(rinv00,_mm_mul_pd(vvdw,dsw)) );
1369             cutoff_mask      = _mm_cmplt_pd(rsq00,rcutoff2);
1370
1371             fscal            = fvdw;
1372
1373             fscal            = _mm_and_pd(fscal,cutoff_mask);
1374
1375             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1376
1377             /* Update vectorial force */
1378             fix0             = _mm_macc_pd(dx00,fscal,fix0);
1379             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
1380             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
1381
1382             fjx0             = _mm_macc_pd(dx00,fscal,fjx0);
1383             fjy0             = _mm_macc_pd(dy00,fscal,fjy0);
1384             fjz0             = _mm_macc_pd(dz00,fscal,fjz0);
1385
1386             }
1387
1388             /**************************
1389              * CALCULATE INTERACTIONS *
1390              **************************/
1391
1392             if (gmx_mm_any_lt(rsq10,rcutoff2))
1393             {
1394
1395             r10              = _mm_mul_pd(rsq10,rinv10);
1396
1397             /* Compute parameters for interactions between i and j atoms */
1398             qq10             = _mm_mul_pd(iq1,jq0);
1399
1400             /* EWALD ELECTROSTATICS */
1401
1402             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1403             ewrt             = _mm_mul_pd(r10,ewtabscale);
1404             ewitab           = _mm_cvttpd_epi32(ewrt);
1405 #ifdef __XOP__
1406             eweps            = _mm_frcz_pd(ewrt);
1407 #else
1408             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1409 #endif
1410             twoeweps         = _mm_add_pd(eweps,eweps);
1411             ewitab           = _mm_slli_epi32(ewitab,2);
1412             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
1413             ewtabD           = _mm_setzero_pd();
1414             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1415             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
1416             ewtabFn          = _mm_setzero_pd();
1417             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1418             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
1419             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
1420             velec            = _mm_mul_pd(qq10,_mm_sub_pd(rinv10,velec));
1421             felec            = _mm_mul_pd(_mm_mul_pd(qq10,rinv10),_mm_sub_pd(rinvsq10,felec));
1422
1423             d                = _mm_sub_pd(r10,rswitch);
1424             d                = _mm_max_pd(d,_mm_setzero_pd());
1425             d2               = _mm_mul_pd(d,d);
1426             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
1427
1428             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
1429
1430             /* Evaluate switch function */
1431             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1432             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv10,_mm_mul_pd(velec,dsw)) );
1433             cutoff_mask      = _mm_cmplt_pd(rsq10,rcutoff2);
1434
1435             fscal            = felec;
1436
1437             fscal            = _mm_and_pd(fscal,cutoff_mask);
1438
1439             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1440
1441             /* Update vectorial force */
1442             fix1             = _mm_macc_pd(dx10,fscal,fix1);
1443             fiy1             = _mm_macc_pd(dy10,fscal,fiy1);
1444             fiz1             = _mm_macc_pd(dz10,fscal,fiz1);
1445
1446             fjx0             = _mm_macc_pd(dx10,fscal,fjx0);
1447             fjy0             = _mm_macc_pd(dy10,fscal,fjy0);
1448             fjz0             = _mm_macc_pd(dz10,fscal,fjz0);
1449
1450             }
1451
1452             /**************************
1453              * CALCULATE INTERACTIONS *
1454              **************************/
1455
1456             if (gmx_mm_any_lt(rsq20,rcutoff2))
1457             {
1458
1459             r20              = _mm_mul_pd(rsq20,rinv20);
1460
1461             /* Compute parameters for interactions between i and j atoms */
1462             qq20             = _mm_mul_pd(iq2,jq0);
1463
1464             /* EWALD ELECTROSTATICS */
1465
1466             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1467             ewrt             = _mm_mul_pd(r20,ewtabscale);
1468             ewitab           = _mm_cvttpd_epi32(ewrt);
1469 #ifdef __XOP__
1470             eweps            = _mm_frcz_pd(ewrt);
1471 #else
1472             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1473 #endif
1474             twoeweps         = _mm_add_pd(eweps,eweps);
1475             ewitab           = _mm_slli_epi32(ewitab,2);
1476             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
1477             ewtabD           = _mm_setzero_pd();
1478             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1479             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
1480             ewtabFn          = _mm_setzero_pd();
1481             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1482             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
1483             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
1484             velec            = _mm_mul_pd(qq20,_mm_sub_pd(rinv20,velec));
1485             felec            = _mm_mul_pd(_mm_mul_pd(qq20,rinv20),_mm_sub_pd(rinvsq20,felec));
1486
1487             d                = _mm_sub_pd(r20,rswitch);
1488             d                = _mm_max_pd(d,_mm_setzero_pd());
1489             d2               = _mm_mul_pd(d,d);
1490             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
1491
1492             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
1493
1494             /* Evaluate switch function */
1495             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1496             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv20,_mm_mul_pd(velec,dsw)) );
1497             cutoff_mask      = _mm_cmplt_pd(rsq20,rcutoff2);
1498
1499             fscal            = felec;
1500
1501             fscal            = _mm_and_pd(fscal,cutoff_mask);
1502
1503             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1504
1505             /* Update vectorial force */
1506             fix2             = _mm_macc_pd(dx20,fscal,fix2);
1507             fiy2             = _mm_macc_pd(dy20,fscal,fiy2);
1508             fiz2             = _mm_macc_pd(dz20,fscal,fiz2);
1509
1510             fjx0             = _mm_macc_pd(dx20,fscal,fjx0);
1511             fjy0             = _mm_macc_pd(dy20,fscal,fjy0);
1512             fjz0             = _mm_macc_pd(dz20,fscal,fjz0);
1513
1514             }
1515
1516             /**************************
1517              * CALCULATE INTERACTIONS *
1518              **************************/
1519
1520             if (gmx_mm_any_lt(rsq30,rcutoff2))
1521             {
1522
1523             r30              = _mm_mul_pd(rsq30,rinv30);
1524
1525             /* Compute parameters for interactions between i and j atoms */
1526             qq30             = _mm_mul_pd(iq3,jq0);
1527
1528             /* EWALD ELECTROSTATICS */
1529
1530             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1531             ewrt             = _mm_mul_pd(r30,ewtabscale);
1532             ewitab           = _mm_cvttpd_epi32(ewrt);
1533 #ifdef __XOP__
1534             eweps            = _mm_frcz_pd(ewrt);
1535 #else
1536             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1537 #endif
1538             twoeweps         = _mm_add_pd(eweps,eweps);
1539             ewitab           = _mm_slli_epi32(ewitab,2);
1540             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
1541             ewtabD           = _mm_setzero_pd();
1542             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1543             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
1544             ewtabFn          = _mm_setzero_pd();
1545             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1546             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
1547             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
1548             velec            = _mm_mul_pd(qq30,_mm_sub_pd(rinv30,velec));
1549             felec            = _mm_mul_pd(_mm_mul_pd(qq30,rinv30),_mm_sub_pd(rinvsq30,felec));
1550
1551             d                = _mm_sub_pd(r30,rswitch);
1552             d                = _mm_max_pd(d,_mm_setzero_pd());
1553             d2               = _mm_mul_pd(d,d);
1554             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
1555
1556             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
1557
1558             /* Evaluate switch function */
1559             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1560             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv30,_mm_mul_pd(velec,dsw)) );
1561             cutoff_mask      = _mm_cmplt_pd(rsq30,rcutoff2);
1562
1563             fscal            = felec;
1564
1565             fscal            = _mm_and_pd(fscal,cutoff_mask);
1566
1567             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1568
1569             /* Update vectorial force */
1570             fix3             = _mm_macc_pd(dx30,fscal,fix3);
1571             fiy3             = _mm_macc_pd(dy30,fscal,fiy3);
1572             fiz3             = _mm_macc_pd(dz30,fscal,fiz3);
1573
1574             fjx0             = _mm_macc_pd(dx30,fscal,fjx0);
1575             fjy0             = _mm_macc_pd(dy30,fscal,fjy0);
1576             fjz0             = _mm_macc_pd(dz30,fscal,fjz0);
1577
1578             }
1579
1580             gmx_mm_decrement_1rvec_1ptr_swizzle_pd(f+j_coord_offsetA,fjx0,fjy0,fjz0);
1581
1582             /* Inner loop uses 257 flops */
1583         }
1584
1585         /* End of innermost loop */
1586
1587         gmx_mm_update_iforce_4atom_swizzle_pd(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
1588                                               f+i_coord_offset,fshift+i_shift_offset);
1589
1590         /* Increment number of inner iterations */
1591         inneriter                  += j_index_end - j_index_start;
1592
1593         /* Outer loop uses 24 flops */
1594     }
1595
1596     /* Increment number of outer iterations */
1597     outeriter        += nri;
1598
1599     /* Update outer/inner flops */
1600
1601     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4_F,outeriter*24 + inneriter*257);
1602 }