src/gromacs/gmxlib/nonbonded/nb_kernel_sse2_double/nb_kernel_ElecRF_VdwCSTab_GeomW4W4_sse2_double.c

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