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