src/gromacs/gmxlib/nonbonded/nb_kernel_sse4_1_double/nb_kernel_ElecCoul_VdwCSTab_GeomW4W4_sse4_1_double.cpp

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