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