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