src/gromacs/gmxlib/nonbonded/nb_kernel_sse4_1_double/nb_kernel_ElecEw_VdwLJ_GeomW3W3_sse4_1_double.c

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