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