src/gromacs/gmxlib/nonbonded/nb_kernel_avx_128_fma_double/nb_kernel_ElecEwSw_VdwLJSw_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 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <math.h>
  43
  44 #include "../nb_kernel.h"
  45 #include "types/simple.h"
  46 #include "gromacs/math/vec.h"
  47 #include "nrnb.h"
  48
  49 #include "gromacs/simd/math_x86_avx_128_fma_double.h"
  50 #include "kernelutil_x86_avx_128_fma_double.h"
  51
  52 /*
  53  * Gromacs nonbonded kernel:   nb_kernel_ElecEwSw_VdwLJSw_GeomW4W4_VF_avx_128_fma_double
  54  * Electrostatics interaction: Ewald
  55  * VdW interaction:            LennardJones
  56  * Geometry:                   Water4-Water4
  57  * Calculate force/pot:        PotentialAndForce
  58  */
  59 void
  60 nb_kernel_ElecEwSw_VdwLJSw_GeomW4W4_VF_avx_128_fma_double
  61                     (t_nblist                    * gmx_restrict       nlist,
  62                      rvec                        * gmx_restrict          xx,
  63                      rvec                        * gmx_restrict          ff,
  64                      t_forcerec                  * gmx_restrict          fr,
  65                      t_mdatoms                   * gmx_restrict     mdatoms,
  66                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
  67                      t_nrnb                      * gmx_restrict        nrnb)
  68 {
  69     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
  70      * just 0 for non-waters.
  71      * Suffixes A,B refer to j loop unrolling done with SSE double precision, e.g. for the two different
  72      * jnr indices corresponding to data put in the four positions in the SIMD register.
  73      */
  74     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
  75     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
  76     int              jnrA,jnrB;
  77     int              j_coord_offsetA,j_coord_offsetB;
  78     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
  79     real             rcutoff_scalar;
  80     real             *shiftvec,*fshift,*x,*f;
  81     __m128d          tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
  82     int              vdwioffset0;
  83     __m128d          ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
  84     int              vdwioffset1;
  85     __m128d          ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
  86     int              vdwioffset2;
  87     __m128d          ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
  88     int              vdwioffset3;
  89     __m128d          ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
  90     int              vdwjidx0A,vdwjidx0B;
  91     __m128d          jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
  92     int              vdwjidx1A,vdwjidx1B;
  93     __m128d          jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
  94     int              vdwjidx2A,vdwjidx2B;
  95     __m128d          jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
  96     int              vdwjidx3A,vdwjidx3B;
  97     __m128d          jx3,jy3,jz3,fjx3,fjy3,fjz3,jq3,isaj3;
  98     __m128d          dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
  99     __m128d          dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11;
 100     __m128d          dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12;
 101     __m128d          dx13,dy13,dz13,rsq13,rinv13,rinvsq13,r13,qq13,c6_13,c12_13;
 102     __m128d          dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21;
 103     __m128d          dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22;
 104     __m128d          dx23,dy23,dz23,rsq23,rinv23,rinvsq23,r23,qq23,c6_23,c12_23;
 105     __m128d          dx31,dy31,dz31,rsq31,rinv31,rinvsq31,r31,qq31,c6_31,c12_31;
 106     __m128d          dx32,dy32,dz32,rsq32,rinv32,rinvsq32,r32,qq32,c6_32,c12_32;
 107     __m128d          dx33,dy33,dz33,rsq33,rinv33,rinvsq33,r33,qq33,c6_33,c12_33;
 108     __m128d          velec,felec,velecsum,facel,crf,krf,krf2;
 109     real             *charge;
 110     int              nvdwtype;
 111     __m128d          rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
 112     int              *vdwtype;
 113     real             *vdwparam;
 114     __m128d          one_sixth   = _mm_set1_pd(1.0/6.0);
 115     __m128d          one_twelfth = _mm_set1_pd(1.0/12.0);
 116     __m128i          ewitab;
 117     __m128d          ewtabscale,eweps,twoeweps,sh_ewald,ewrt,ewtabhalfspace,ewtabF,ewtabFn,ewtabD,ewtabV;
 118     real             *ewtab;
 119     __m128d          rswitch,swV3,swV4,swV5,swF2,swF3,swF4,d,d2,sw,dsw;
 120     real             rswitch_scalar,d_scalar;
 121     __m128d          dummy_mask,cutoff_mask;
 122     __m128d          signbit   = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
 123     __m128d          one     = _mm_set1_pd(1.0);
 124     __m128d          two     = _mm_set1_pd(2.0);
 125     x                = xx[0];
 126     f                = ff[0];
 127
 128     nri              = nlist->nri;
 129     iinr             = nlist->iinr;
 130     jindex           = nlist->jindex;
 131     jjnr             = nlist->jjnr;
 132     shiftidx         = nlist->shift;
 133     gid              = nlist->gid;
 134     shiftvec         = fr->shift_vec[0];
 135     fshift           = fr->fshift[0];
 136     facel            = _mm_set1_pd(fr->epsfac);
 137     charge           = mdatoms->chargeA;
 138     nvdwtype         = fr->ntype;
 139     vdwparam         = fr->nbfp;
 140     vdwtype          = mdatoms->typeA;
 141
 142     sh_ewald         = _mm_set1_pd(fr->ic->sh_ewald);
 143     ewtab            = fr->ic->tabq_coul_FDV0;
 144     ewtabscale       = _mm_set1_pd(fr->ic->tabq_scale);
 145     ewtabhalfspace   = _mm_set1_pd(0.5/fr->ic->tabq_scale);
 146
 147     /* Setup water-specific parameters */
 148     inr              = nlist->iinr[0];
 149     iq1              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+1]));
 150     iq2              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+2]));
 151     iq3              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+3]));
 152     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
 153
 154     jq1              = _mm_set1_pd(charge[inr+1]);
 155     jq2              = _mm_set1_pd(charge[inr+2]);
 156     jq3              = _mm_set1_pd(charge[inr+3]);
 157     vdwjidx0A        = 2*vdwtype[inr+0];
 158     c6_00            = _mm_set1_pd(vdwparam[vdwioffset0+vdwjidx0A]);
 159     c12_00           = _mm_set1_pd(vdwparam[vdwioffset0+vdwjidx0A+1]);
 160     qq11             = _mm_mul_pd(iq1,jq1);
 161     qq12             = _mm_mul_pd(iq1,jq2);
 162     qq13             = _mm_mul_pd(iq1,jq3);
 163     qq21             = _mm_mul_pd(iq2,jq1);
 164     qq22             = _mm_mul_pd(iq2,jq2);
 165     qq23             = _mm_mul_pd(iq2,jq3);
 166     qq31             = _mm_mul_pd(iq3,jq1);
 167     qq32             = _mm_mul_pd(iq3,jq2);
 168     qq33             = _mm_mul_pd(iq3,jq3);
 169
 170     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
 171     rcutoff_scalar   = fr->rcoulomb;
 172     rcutoff          = _mm_set1_pd(rcutoff_scalar);
 173     rcutoff2         = _mm_mul_pd(rcutoff,rcutoff);
 174
 175     rswitch_scalar   = fr->rcoulomb_switch;
 176     rswitch          = _mm_set1_pd(rswitch_scalar);
 177     /* Setup switch parameters */
 178     d_scalar         = rcutoff_scalar-rswitch_scalar;
 179     d                = _mm_set1_pd(d_scalar);
 180     swV3             = _mm_set1_pd(-10.0/(d_scalar*d_scalar*d_scalar));
 181     swV4             = _mm_set1_pd( 15.0/(d_scalar*d_scalar*d_scalar*d_scalar));
 182     swV5             = _mm_set1_pd( -6.0/(d_scalar*d_scalar*d_scalar*d_scalar*d_scalar));
 183     swF2             = _mm_set1_pd(-30.0/(d_scalar*d_scalar*d_scalar));
 184     swF3             = _mm_set1_pd( 60.0/(d_scalar*d_scalar*d_scalar*d_scalar));
 185     swF4             = _mm_set1_pd(-30.0/(d_scalar*d_scalar*d_scalar*d_scalar*d_scalar));
 186
 187     /* Avoid stupid compiler warnings */
 188     jnrA = jnrB = 0;
 189     j_coord_offsetA = 0;
 190     j_coord_offsetB = 0;
 191
 192     outeriter        = 0;
 193     inneriter        = 0;
 194
 195     /* Start outer loop over neighborlists */
 196     for(iidx=0; iidx<nri; iidx++)
 197     {
 198         /* Load shift vector for this list */
 199         i_shift_offset   = DIM*shiftidx[iidx];
 200
 201         /* Load limits for loop over neighbors */
 202         j_index_start    = jindex[iidx];
 203         j_index_end      = jindex[iidx+1];
 204
 205         /* Get outer coordinate index */
 206         inr              = iinr[iidx];
 207         i_coord_offset   = DIM*inr;
 208
 209         /* Load i particle coords and add shift vector */
 210         gmx_mm_load_shift_and_4rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,
 211                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
 212
 213         fix0             = _mm_setzero_pd();
 214         fiy0             = _mm_setzero_pd();
 215         fiz0             = _mm_setzero_pd();
 216         fix1             = _mm_setzero_pd();
 217         fiy1             = _mm_setzero_pd();
 218         fiz1             = _mm_setzero_pd();
 219         fix2             = _mm_setzero_pd();
 220         fiy2             = _mm_setzero_pd();
 221         fiz2             = _mm_setzero_pd();
 222         fix3             = _mm_setzero_pd();
 223         fiy3             = _mm_setzero_pd();
 224         fiz3             = _mm_setzero_pd();
 225
 226         /* Reset potential sums */
 227         velecsum         = _mm_setzero_pd();
 228         vvdwsum          = _mm_setzero_pd();
 229
 230         /* Start inner kernel loop */
 231         for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
 232         {
 233
 234             /* Get j neighbor index, and coordinate index */
 235             jnrA             = jjnr[jidx];
 236             jnrB             = jjnr[jidx+1];
 237             j_coord_offsetA  = DIM*jnrA;
 238             j_coord_offsetB  = DIM*jnrB;
 239
 240             /* load j atom coordinates */
 241             gmx_mm_load_4rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
 242                                               &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,
 243                                               &jy2,&jz2,&jx3,&jy3,&jz3);
 244
 245             /* Calculate displacement vector */
 246             dx00             = _mm_sub_pd(ix0,jx0);
 247             dy00             = _mm_sub_pd(iy0,jy0);
 248             dz00             = _mm_sub_pd(iz0,jz0);
 249             dx11             = _mm_sub_pd(ix1,jx1);
 250             dy11             = _mm_sub_pd(iy1,jy1);
 251             dz11             = _mm_sub_pd(iz1,jz1);
 252             dx12             = _mm_sub_pd(ix1,jx2);
 253             dy12             = _mm_sub_pd(iy1,jy2);
 254             dz12             = _mm_sub_pd(iz1,jz2);
 255             dx13             = _mm_sub_pd(ix1,jx3);
 256             dy13             = _mm_sub_pd(iy1,jy3);
 257             dz13             = _mm_sub_pd(iz1,jz3);
 258             dx21             = _mm_sub_pd(ix2,jx1);
 259             dy21             = _mm_sub_pd(iy2,jy1);
 260             dz21             = _mm_sub_pd(iz2,jz1);
 261             dx22             = _mm_sub_pd(ix2,jx2);
 262             dy22             = _mm_sub_pd(iy2,jy2);
 263             dz22             = _mm_sub_pd(iz2,jz2);
 264             dx23             = _mm_sub_pd(ix2,jx3);
 265             dy23             = _mm_sub_pd(iy2,jy3);
 266             dz23             = _mm_sub_pd(iz2,jz3);
 267             dx31             = _mm_sub_pd(ix3,jx1);
 268             dy31             = _mm_sub_pd(iy3,jy1);
 269             dz31             = _mm_sub_pd(iz3,jz1);
 270             dx32             = _mm_sub_pd(ix3,jx2);
 271             dy32             = _mm_sub_pd(iy3,jy2);
 272             dz32             = _mm_sub_pd(iz3,jz2);
 273             dx33             = _mm_sub_pd(ix3,jx3);
 274             dy33             = _mm_sub_pd(iy3,jy3);
 275             dz33             = _mm_sub_pd(iz3,jz3);
 276
 277             /* Calculate squared distance and things based on it */
 278             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
 279             rsq11            = gmx_mm_calc_rsq_pd(dx11,dy11,dz11);
 280             rsq12            = gmx_mm_calc_rsq_pd(dx12,dy12,dz12);
 281             rsq13            = gmx_mm_calc_rsq_pd(dx13,dy13,dz13);
 282             rsq21            = gmx_mm_calc_rsq_pd(dx21,dy21,dz21);
 283             rsq22            = gmx_mm_calc_rsq_pd(dx22,dy22,dz22);
 284             rsq23            = gmx_mm_calc_rsq_pd(dx23,dy23,dz23);
 285             rsq31            = gmx_mm_calc_rsq_pd(dx31,dy31,dz31);
 286             rsq32            = gmx_mm_calc_rsq_pd(dx32,dy32,dz32);
 287             rsq33            = gmx_mm_calc_rsq_pd(dx33,dy33,dz33);
 288
 289             rinv00           = gmx_mm_invsqrt_pd(rsq00);
 290             rinv11           = gmx_mm_invsqrt_pd(rsq11);
 291             rinv12           = gmx_mm_invsqrt_pd(rsq12);
 292             rinv13           = gmx_mm_invsqrt_pd(rsq13);
 293             rinv21           = gmx_mm_invsqrt_pd(rsq21);
 294             rinv22           = gmx_mm_invsqrt_pd(rsq22);
 295             rinv23           = gmx_mm_invsqrt_pd(rsq23);
 296             rinv31           = gmx_mm_invsqrt_pd(rsq31);
 297             rinv32           = gmx_mm_invsqrt_pd(rsq32);
 298             rinv33           = gmx_mm_invsqrt_pd(rsq33);
 299
 300             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
 301             rinvsq11         = _mm_mul_pd(rinv11,rinv11);
 302             rinvsq12         = _mm_mul_pd(rinv12,rinv12);
 303             rinvsq13         = _mm_mul_pd(rinv13,rinv13);
 304             rinvsq21         = _mm_mul_pd(rinv21,rinv21);
 305             rinvsq22         = _mm_mul_pd(rinv22,rinv22);
 306             rinvsq23         = _mm_mul_pd(rinv23,rinv23);
 307             rinvsq31         = _mm_mul_pd(rinv31,rinv31);
 308             rinvsq32         = _mm_mul_pd(rinv32,rinv32);
 309             rinvsq33         = _mm_mul_pd(rinv33,rinv33);
 310
 311             fjx0             = _mm_setzero_pd();
 312             fjy0             = _mm_setzero_pd();
 313             fjz0             = _mm_setzero_pd();
 314             fjx1             = _mm_setzero_pd();
 315             fjy1             = _mm_setzero_pd();
 316             fjz1             = _mm_setzero_pd();
 317             fjx2             = _mm_setzero_pd();
 318             fjy2             = _mm_setzero_pd();
 319             fjz2             = _mm_setzero_pd();
 320             fjx3             = _mm_setzero_pd();
 321             fjy3             = _mm_setzero_pd();
 322             fjz3             = _mm_setzero_pd();
 323
 324             /**************************
 325              * CALCULATE INTERACTIONS *
 326              **************************/
 327
 328             if (gmx_mm_any_lt(rsq00,rcutoff2))
 329             {
 330
 331             r00              = _mm_mul_pd(rsq00,rinv00);
 332
 333             /* LENNARD-JONES DISPERSION/REPULSION */
 334
 335             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
 336             vvdw6            = _mm_mul_pd(c6_00,rinvsix);
 337             vvdw12           = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
 338             vvdw             = _mm_msub_pd( vvdw12,one_twelfth, _mm_mul_pd(vvdw6,one_sixth) );
 339             fvdw             = _mm_mul_pd(_mm_sub_pd(vvdw12,vvdw6),rinvsq00);
 340
 341             d                = _mm_sub_pd(r00,rswitch);
 342             d                = _mm_max_pd(d,_mm_setzero_pd());
 343             d2               = _mm_mul_pd(d,d);
 344             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
 345
 346             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
 347
 348             /* Evaluate switch function */
 349             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 350             fvdw             = _mm_msub_pd( fvdw,sw , _mm_mul_pd(rinv00,_mm_mul_pd(vvdw,dsw)) );
 351             vvdw             = _mm_mul_pd(vvdw,sw);
 352             cutoff_mask      = _mm_cmplt_pd(rsq00,rcutoff2);
 353
 354             /* Update potential sum for this i atom from the interaction with this j atom. */
 355             vvdw             = _mm_and_pd(vvdw,cutoff_mask);
 356             vvdwsum          = _mm_add_pd(vvdwsum,vvdw);
 357
 358             fscal            = fvdw;
 359
 360             fscal            = _mm_and_pd(fscal,cutoff_mask);
 361
 362             /* Update vectorial force */
 363             fix0             = _mm_macc_pd(dx00,fscal,fix0);
 364             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
 365             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
 366
 367             fjx0             = _mm_macc_pd(dx00,fscal,fjx0);
 368             fjy0             = _mm_macc_pd(dy00,fscal,fjy0);
 369             fjz0             = _mm_macc_pd(dz00,fscal,fjz0);
 370
 371             }
 372
 373             /**************************
 374              * CALCULATE INTERACTIONS *
 375              **************************/
 376
 377             if (gmx_mm_any_lt(rsq11,rcutoff2))
 378             {
 379
 380             r11              = _mm_mul_pd(rsq11,rinv11);
 381
 382             /* EWALD ELECTROSTATICS */
 383
 384             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 385             ewrt             = _mm_mul_pd(r11,ewtabscale);
 386             ewitab           = _mm_cvttpd_epi32(ewrt);
 387 #ifdef __XOP__
 388             eweps            = _mm_frcz_pd(ewrt);
 389 #else
 390             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 391 #endif
 392             twoeweps         = _mm_add_pd(eweps,eweps);
 393             ewitab           = _mm_slli_epi32(ewitab,2);
 394             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 395             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 396             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 397             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 398             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
 399             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 400             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 401             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 402             velec            = _mm_mul_pd(qq11,_mm_sub_pd(rinv11,velec));
 403             felec            = _mm_mul_pd(_mm_mul_pd(qq11,rinv11),_mm_sub_pd(rinvsq11,felec));
 404
 405             d                = _mm_sub_pd(r11,rswitch);
 406             d                = _mm_max_pd(d,_mm_setzero_pd());
 407             d2               = _mm_mul_pd(d,d);
 408             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
 409
 410             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
 411
 412             /* Evaluate switch function */
 413             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 414             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv11,_mm_mul_pd(velec,dsw)) );
 415             velec            = _mm_mul_pd(velec,sw);
 416             cutoff_mask      = _mm_cmplt_pd(rsq11,rcutoff2);
 417
 418             /* Update potential sum for this i atom from the interaction with this j atom. */
 419             velec            = _mm_and_pd(velec,cutoff_mask);
 420             velecsum         = _mm_add_pd(velecsum,velec);
 421
 422             fscal            = felec;
 423
 424             fscal            = _mm_and_pd(fscal,cutoff_mask);
 425
 426             /* Update vectorial force */
 427             fix1             = _mm_macc_pd(dx11,fscal,fix1);
 428             fiy1             = _mm_macc_pd(dy11,fscal,fiy1);
 429             fiz1             = _mm_macc_pd(dz11,fscal,fiz1);
 430
 431             fjx1             = _mm_macc_pd(dx11,fscal,fjx1);
 432             fjy1             = _mm_macc_pd(dy11,fscal,fjy1);
 433             fjz1             = _mm_macc_pd(dz11,fscal,fjz1);
 434
 435             }
 436
 437             /**************************
 438              * CALCULATE INTERACTIONS *
 439              **************************/
 440
 441             if (gmx_mm_any_lt(rsq12,rcutoff2))
 442             {
 443
 444             r12              = _mm_mul_pd(rsq12,rinv12);
 445
 446             /* EWALD ELECTROSTATICS */
 447
 448             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 449             ewrt             = _mm_mul_pd(r12,ewtabscale);
 450             ewitab           = _mm_cvttpd_epi32(ewrt);
 451 #ifdef __XOP__
 452             eweps            = _mm_frcz_pd(ewrt);
 453 #else
 454             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 455 #endif
 456             twoeweps         = _mm_add_pd(eweps,eweps);
 457             ewitab           = _mm_slli_epi32(ewitab,2);
 458             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 459             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 460             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 461             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 462             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
 463             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 464             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 465             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 466             velec            = _mm_mul_pd(qq12,_mm_sub_pd(rinv12,velec));
 467             felec            = _mm_mul_pd(_mm_mul_pd(qq12,rinv12),_mm_sub_pd(rinvsq12,felec));
 468
 469             d                = _mm_sub_pd(r12,rswitch);
 470             d                = _mm_max_pd(d,_mm_setzero_pd());
 471             d2               = _mm_mul_pd(d,d);
 472             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
 473
 474             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
 475
 476             /* Evaluate switch function */
 477             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 478             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv12,_mm_mul_pd(velec,dsw)) );
 479             velec            = _mm_mul_pd(velec,sw);
 480             cutoff_mask      = _mm_cmplt_pd(rsq12,rcutoff2);
 481
 482             /* Update potential sum for this i atom from the interaction with this j atom. */
 483             velec            = _mm_and_pd(velec,cutoff_mask);
 484             velecsum         = _mm_add_pd(velecsum,velec);
 485
 486             fscal            = felec;
 487
 488             fscal            = _mm_and_pd(fscal,cutoff_mask);
 489
 490             /* Update vectorial force */
 491             fix1             = _mm_macc_pd(dx12,fscal,fix1);
 492             fiy1             = _mm_macc_pd(dy12,fscal,fiy1);
 493             fiz1             = _mm_macc_pd(dz12,fscal,fiz1);
 494
 495             fjx2             = _mm_macc_pd(dx12,fscal,fjx2);
 496             fjy2             = _mm_macc_pd(dy12,fscal,fjy2);
 497             fjz2             = _mm_macc_pd(dz12,fscal,fjz2);
 498
 499             }
 500
 501             /**************************
 502              * CALCULATE INTERACTIONS *
 503              **************************/
 504
 505             if (gmx_mm_any_lt(rsq13,rcutoff2))
 506             {
 507
 508             r13              = _mm_mul_pd(rsq13,rinv13);
 509
 510             /* EWALD ELECTROSTATICS */
 511
 512             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 513             ewrt             = _mm_mul_pd(r13,ewtabscale);
 514             ewitab           = _mm_cvttpd_epi32(ewrt);
 515 #ifdef __XOP__
 516             eweps            = _mm_frcz_pd(ewrt);
 517 #else
 518             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 519 #endif
 520             twoeweps         = _mm_add_pd(eweps,eweps);
 521             ewitab           = _mm_slli_epi32(ewitab,2);
 522             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 523             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 524             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 525             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 526             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
 527             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 528             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 529             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 530             velec            = _mm_mul_pd(qq13,_mm_sub_pd(rinv13,velec));
 531             felec            = _mm_mul_pd(_mm_mul_pd(qq13,rinv13),_mm_sub_pd(rinvsq13,felec));
 532
 533             d                = _mm_sub_pd(r13,rswitch);
 534             d                = _mm_max_pd(d,_mm_setzero_pd());
 535             d2               = _mm_mul_pd(d,d);
 536             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
 537
 538             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
 539
 540             /* Evaluate switch function */
 541             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 542             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv13,_mm_mul_pd(velec,dsw)) );
 543             velec            = _mm_mul_pd(velec,sw);
 544             cutoff_mask      = _mm_cmplt_pd(rsq13,rcutoff2);
 545
 546             /* Update potential sum for this i atom from the interaction with this j atom. */
 547             velec            = _mm_and_pd(velec,cutoff_mask);
 548             velecsum         = _mm_add_pd(velecsum,velec);
 549
 550             fscal            = felec;
 551
 552             fscal            = _mm_and_pd(fscal,cutoff_mask);
 553
 554             /* Update vectorial force */
 555             fix1             = _mm_macc_pd(dx13,fscal,fix1);
 556             fiy1             = _mm_macc_pd(dy13,fscal,fiy1);
 557             fiz1             = _mm_macc_pd(dz13,fscal,fiz1);
 558
 559             fjx3             = _mm_macc_pd(dx13,fscal,fjx3);
 560             fjy3             = _mm_macc_pd(dy13,fscal,fjy3);
 561             fjz3             = _mm_macc_pd(dz13,fscal,fjz3);
 562
 563             }
 564
 565             /**************************
 566              * CALCULATE INTERACTIONS *
 567              **************************/
 568
 569             if (gmx_mm_any_lt(rsq21,rcutoff2))
 570             {
 571
 572             r21              = _mm_mul_pd(rsq21,rinv21);
 573
 574             /* EWALD ELECTROSTATICS */
 575
 576             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 577             ewrt             = _mm_mul_pd(r21,ewtabscale);
 578             ewitab           = _mm_cvttpd_epi32(ewrt);
 579 #ifdef __XOP__
 580             eweps            = _mm_frcz_pd(ewrt);
 581 #else
 582             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 583 #endif
 584             twoeweps         = _mm_add_pd(eweps,eweps);
 585             ewitab           = _mm_slli_epi32(ewitab,2);
 586             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 587             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 588             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 589             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 590             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
 591             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 592             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 593             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 594             velec            = _mm_mul_pd(qq21,_mm_sub_pd(rinv21,velec));
 595             felec            = _mm_mul_pd(_mm_mul_pd(qq21,rinv21),_mm_sub_pd(rinvsq21,felec));
 596
 597             d                = _mm_sub_pd(r21,rswitch);
 598             d                = _mm_max_pd(d,_mm_setzero_pd());
 599             d2               = _mm_mul_pd(d,d);
 600             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
 601
 602             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
 603
 604             /* Evaluate switch function */
 605             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 606             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv21,_mm_mul_pd(velec,dsw)) );
 607             velec            = _mm_mul_pd(velec,sw);
 608             cutoff_mask      = _mm_cmplt_pd(rsq21,rcutoff2);
 609
 610             /* Update potential sum for this i atom from the interaction with this j atom. */
 611             velec            = _mm_and_pd(velec,cutoff_mask);
 612             velecsum         = _mm_add_pd(velecsum,velec);
 613
 614             fscal            = felec;
 615
 616             fscal            = _mm_and_pd(fscal,cutoff_mask);
 617
 618             /* Update vectorial force */
 619             fix2             = _mm_macc_pd(dx21,fscal,fix2);
 620             fiy2             = _mm_macc_pd(dy21,fscal,fiy2);
 621             fiz2             = _mm_macc_pd(dz21,fscal,fiz2);
 622
 623             fjx1             = _mm_macc_pd(dx21,fscal,fjx1);
 624             fjy1             = _mm_macc_pd(dy21,fscal,fjy1);
 625             fjz1             = _mm_macc_pd(dz21,fscal,fjz1);
 626
 627             }
 628
 629             /**************************
 630              * CALCULATE INTERACTIONS *
 631              **************************/
 632
 633             if (gmx_mm_any_lt(rsq22,rcutoff2))
 634             {
 635
 636             r22              = _mm_mul_pd(rsq22,rinv22);
 637
 638             /* EWALD ELECTROSTATICS */
 639
 640             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 641             ewrt             = _mm_mul_pd(r22,ewtabscale);
 642             ewitab           = _mm_cvttpd_epi32(ewrt);
 643 #ifdef __XOP__
 644             eweps            = _mm_frcz_pd(ewrt);
 645 #else
 646             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 647 #endif
 648             twoeweps         = _mm_add_pd(eweps,eweps);
 649             ewitab           = _mm_slli_epi32(ewitab,2);
 650             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 651             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 652             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 653             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 654             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
 655             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 656             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 657             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 658             velec            = _mm_mul_pd(qq22,_mm_sub_pd(rinv22,velec));
 659             felec            = _mm_mul_pd(_mm_mul_pd(qq22,rinv22),_mm_sub_pd(rinvsq22,felec));
 660
 661             d                = _mm_sub_pd(r22,rswitch);
 662             d                = _mm_max_pd(d,_mm_setzero_pd());
 663             d2               = _mm_mul_pd(d,d);
 664             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
 665
 666             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
 667
 668             /* Evaluate switch function */
 669             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 670             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv22,_mm_mul_pd(velec,dsw)) );
 671             velec            = _mm_mul_pd(velec,sw);
 672             cutoff_mask      = _mm_cmplt_pd(rsq22,rcutoff2);
 673
 674             /* Update potential sum for this i atom from the interaction with this j atom. */
 675             velec            = _mm_and_pd(velec,cutoff_mask);
 676             velecsum         = _mm_add_pd(velecsum,velec);
 677
 678             fscal            = felec;
 679
 680             fscal            = _mm_and_pd(fscal,cutoff_mask);
 681
 682             /* Update vectorial force */
 683             fix2             = _mm_macc_pd(dx22,fscal,fix2);
 684             fiy2             = _mm_macc_pd(dy22,fscal,fiy2);
 685             fiz2             = _mm_macc_pd(dz22,fscal,fiz2);
 686
 687             fjx2             = _mm_macc_pd(dx22,fscal,fjx2);
 688             fjy2             = _mm_macc_pd(dy22,fscal,fjy2);
 689             fjz2             = _mm_macc_pd(dz22,fscal,fjz2);
 690
 691             }
 692
 693             /**************************
 694              * CALCULATE INTERACTIONS *
 695              **************************/
 696
 697             if (gmx_mm_any_lt(rsq23,rcutoff2))
 698             {
 699
 700             r23              = _mm_mul_pd(rsq23,rinv23);
 701
 702             /* EWALD ELECTROSTATICS */
 703
 704             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 705             ewrt             = _mm_mul_pd(r23,ewtabscale);
 706             ewitab           = _mm_cvttpd_epi32(ewrt);
 707 #ifdef __XOP__
 708             eweps            = _mm_frcz_pd(ewrt);
 709 #else
 710             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 711 #endif
 712             twoeweps         = _mm_add_pd(eweps,eweps);
 713             ewitab           = _mm_slli_epi32(ewitab,2);
 714             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 715             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 716             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 717             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 718             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
 719             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 720             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 721             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 722             velec            = _mm_mul_pd(qq23,_mm_sub_pd(rinv23,velec));
 723             felec            = _mm_mul_pd(_mm_mul_pd(qq23,rinv23),_mm_sub_pd(rinvsq23,felec));
 724
 725             d                = _mm_sub_pd(r23,rswitch);
 726             d                = _mm_max_pd(d,_mm_setzero_pd());
 727             d2               = _mm_mul_pd(d,d);
 728             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
 729
 730             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
 731
 732             /* Evaluate switch function */
 733             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 734             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv23,_mm_mul_pd(velec,dsw)) );
 735             velec            = _mm_mul_pd(velec,sw);
 736             cutoff_mask      = _mm_cmplt_pd(rsq23,rcutoff2);
 737
 738             /* Update potential sum for this i atom from the interaction with this j atom. */
 739             velec            = _mm_and_pd(velec,cutoff_mask);
 740             velecsum         = _mm_add_pd(velecsum,velec);
 741
 742             fscal            = felec;
 743
 744             fscal            = _mm_and_pd(fscal,cutoff_mask);
 745
 746             /* Update vectorial force */
 747             fix2             = _mm_macc_pd(dx23,fscal,fix2);
 748             fiy2             = _mm_macc_pd(dy23,fscal,fiy2);
 749             fiz2             = _mm_macc_pd(dz23,fscal,fiz2);
 750
 751             fjx3             = _mm_macc_pd(dx23,fscal,fjx3);
 752             fjy3             = _mm_macc_pd(dy23,fscal,fjy3);
 753             fjz3             = _mm_macc_pd(dz23,fscal,fjz3);
 754
 755             }
 756
 757             /**************************
 758              * CALCULATE INTERACTIONS *
 759              **************************/
 760
 761             if (gmx_mm_any_lt(rsq31,rcutoff2))
 762             {
 763
 764             r31              = _mm_mul_pd(rsq31,rinv31);
 765
 766             /* EWALD ELECTROSTATICS */
 767
 768             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 769             ewrt             = _mm_mul_pd(r31,ewtabscale);
 770             ewitab           = _mm_cvttpd_epi32(ewrt);
 771 #ifdef __XOP__
 772             eweps            = _mm_frcz_pd(ewrt);
 773 #else
 774             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 775 #endif
 776             twoeweps         = _mm_add_pd(eweps,eweps);
 777             ewitab           = _mm_slli_epi32(ewitab,2);
 778             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 779             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 780             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 781             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 782             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
 783             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 784             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 785             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 786             velec            = _mm_mul_pd(qq31,_mm_sub_pd(rinv31,velec));
 787             felec            = _mm_mul_pd(_mm_mul_pd(qq31,rinv31),_mm_sub_pd(rinvsq31,felec));
 788
 789             d                = _mm_sub_pd(r31,rswitch);
 790             d                = _mm_max_pd(d,_mm_setzero_pd());
 791             d2               = _mm_mul_pd(d,d);
 792             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
 793
 794             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
 795
 796             /* Evaluate switch function */
 797             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 798             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv31,_mm_mul_pd(velec,dsw)) );
 799             velec            = _mm_mul_pd(velec,sw);
 800             cutoff_mask      = _mm_cmplt_pd(rsq31,rcutoff2);
 801
 802             /* Update potential sum for this i atom from the interaction with this j atom. */
 803             velec            = _mm_and_pd(velec,cutoff_mask);
 804             velecsum         = _mm_add_pd(velecsum,velec);
 805
 806             fscal            = felec;
 807
 808             fscal            = _mm_and_pd(fscal,cutoff_mask);
 809
 810             /* Update vectorial force */
 811             fix3             = _mm_macc_pd(dx31,fscal,fix3);
 812             fiy3             = _mm_macc_pd(dy31,fscal,fiy3);
 813             fiz3             = _mm_macc_pd(dz31,fscal,fiz3);
 814
 815             fjx1             = _mm_macc_pd(dx31,fscal,fjx1);
 816             fjy1             = _mm_macc_pd(dy31,fscal,fjy1);
 817             fjz1             = _mm_macc_pd(dz31,fscal,fjz1);
 818
 819             }
 820
 821             /**************************
 822              * CALCULATE INTERACTIONS *
 823              **************************/
 824
 825             if (gmx_mm_any_lt(rsq32,rcutoff2))
 826             {
 827
 828             r32              = _mm_mul_pd(rsq32,rinv32);
 829
 830             /* EWALD ELECTROSTATICS */
 831
 832             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 833             ewrt             = _mm_mul_pd(r32,ewtabscale);
 834             ewitab           = _mm_cvttpd_epi32(ewrt);
 835 #ifdef __XOP__
 836             eweps            = _mm_frcz_pd(ewrt);
 837 #else
 838             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 839 #endif
 840             twoeweps         = _mm_add_pd(eweps,eweps);
 841             ewitab           = _mm_slli_epi32(ewitab,2);
 842             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 843             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 844             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 845             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 846             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
 847             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 848             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 849             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 850             velec            = _mm_mul_pd(qq32,_mm_sub_pd(rinv32,velec));
 851             felec            = _mm_mul_pd(_mm_mul_pd(qq32,rinv32),_mm_sub_pd(rinvsq32,felec));
 852
 853             d                = _mm_sub_pd(r32,rswitch);
 854             d                = _mm_max_pd(d,_mm_setzero_pd());
 855             d2               = _mm_mul_pd(d,d);
 856             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
 857
 858             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
 859
 860             /* Evaluate switch function */
 861             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 862             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv32,_mm_mul_pd(velec,dsw)) );
 863             velec            = _mm_mul_pd(velec,sw);
 864             cutoff_mask      = _mm_cmplt_pd(rsq32,rcutoff2);
 865
 866             /* Update potential sum for this i atom from the interaction with this j atom. */
 867             velec            = _mm_and_pd(velec,cutoff_mask);
 868             velecsum         = _mm_add_pd(velecsum,velec);
 869
 870             fscal            = felec;
 871
 872             fscal            = _mm_and_pd(fscal,cutoff_mask);
 873
 874             /* Update vectorial force */
 875             fix3             = _mm_macc_pd(dx32,fscal,fix3);
 876             fiy3             = _mm_macc_pd(dy32,fscal,fiy3);
 877             fiz3             = _mm_macc_pd(dz32,fscal,fiz3);
 878
 879             fjx2             = _mm_macc_pd(dx32,fscal,fjx2);
 880             fjy2             = _mm_macc_pd(dy32,fscal,fjy2);
 881             fjz2             = _mm_macc_pd(dz32,fscal,fjz2);
 882
 883             }
 884
 885             /**************************
 886              * CALCULATE INTERACTIONS *
 887              **************************/
 888
 889             if (gmx_mm_any_lt(rsq33,rcutoff2))
 890             {
 891
 892             r33              = _mm_mul_pd(rsq33,rinv33);
 893
 894             /* EWALD ELECTROSTATICS */
 895
 896             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 897             ewrt             = _mm_mul_pd(r33,ewtabscale);
 898             ewitab           = _mm_cvttpd_epi32(ewrt);
 899 #ifdef __XOP__
 900             eweps            = _mm_frcz_pd(ewrt);
 901 #else
 902             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 903 #endif
 904             twoeweps         = _mm_add_pd(eweps,eweps);
 905             ewitab           = _mm_slli_epi32(ewitab,2);
 906             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 907             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 908             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 909             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 910             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
 911             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 912             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 913             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 914             velec            = _mm_mul_pd(qq33,_mm_sub_pd(rinv33,velec));
 915             felec            = _mm_mul_pd(_mm_mul_pd(qq33,rinv33),_mm_sub_pd(rinvsq33,felec));
 916
 917             d                = _mm_sub_pd(r33,rswitch);
 918             d                = _mm_max_pd(d,_mm_setzero_pd());
 919             d2               = _mm_mul_pd(d,d);
 920             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
 921
 922             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
 923
 924             /* Evaluate switch function */
 925             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 926             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv33,_mm_mul_pd(velec,dsw)) );
 927             velec            = _mm_mul_pd(velec,sw);
 928             cutoff_mask      = _mm_cmplt_pd(rsq33,rcutoff2);
 929
 930             /* Update potential sum for this i atom from the interaction with this j atom. */
 931             velec            = _mm_and_pd(velec,cutoff_mask);
 932             velecsum         = _mm_add_pd(velecsum,velec);
 933
 934             fscal            = felec;
 935
 936             fscal            = _mm_and_pd(fscal,cutoff_mask);
 937
 938             /* Update vectorial force */
 939             fix3             = _mm_macc_pd(dx33,fscal,fix3);
 940             fiy3             = _mm_macc_pd(dy33,fscal,fiy3);
 941             fiz3             = _mm_macc_pd(dz33,fscal,fiz3);
 942
 943             fjx3             = _mm_macc_pd(dx33,fscal,fjx3);
 944             fjy3             = _mm_macc_pd(dy33,fscal,fjy3);
 945             fjz3             = _mm_macc_pd(dz33,fscal,fjz3);
 946
 947             }
 948
 949             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);
 950
 951             /* Inner loop uses 677 flops */
 952         }
 953
 954         if(jidx<j_index_end)
 955         {
 956
 957             jnrA             = jjnr[jidx];
 958             j_coord_offsetA  = DIM*jnrA;
 959
 960             /* load j atom coordinates */
 961             gmx_mm_load_4rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
 962                                               &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,
 963                                               &jy2,&jz2,&jx3,&jy3,&jz3);
 964
 965             /* Calculate displacement vector */
 966             dx00             = _mm_sub_pd(ix0,jx0);
 967             dy00             = _mm_sub_pd(iy0,jy0);
 968             dz00             = _mm_sub_pd(iz0,jz0);
 969             dx11             = _mm_sub_pd(ix1,jx1);
 970             dy11             = _mm_sub_pd(iy1,jy1);
 971             dz11             = _mm_sub_pd(iz1,jz1);
 972             dx12             = _mm_sub_pd(ix1,jx2);
 973             dy12             = _mm_sub_pd(iy1,jy2);
 974             dz12             = _mm_sub_pd(iz1,jz2);
 975             dx13             = _mm_sub_pd(ix1,jx3);
 976             dy13             = _mm_sub_pd(iy1,jy3);
 977             dz13             = _mm_sub_pd(iz1,jz3);
 978             dx21             = _mm_sub_pd(ix2,jx1);
 979             dy21             = _mm_sub_pd(iy2,jy1);
 980             dz21             = _mm_sub_pd(iz2,jz1);
 981             dx22             = _mm_sub_pd(ix2,jx2);
 982             dy22             = _mm_sub_pd(iy2,jy2);
 983             dz22             = _mm_sub_pd(iz2,jz2);
 984             dx23             = _mm_sub_pd(ix2,jx3);
 985             dy23             = _mm_sub_pd(iy2,jy3);
 986             dz23             = _mm_sub_pd(iz2,jz3);
 987             dx31             = _mm_sub_pd(ix3,jx1);
 988             dy31             = _mm_sub_pd(iy3,jy1);
 989             dz31             = _mm_sub_pd(iz3,jz1);
 990             dx32             = _mm_sub_pd(ix3,jx2);
 991             dy32             = _mm_sub_pd(iy3,jy2);
 992             dz32             = _mm_sub_pd(iz3,jz2);
 993             dx33             = _mm_sub_pd(ix3,jx3);
 994             dy33             = _mm_sub_pd(iy3,jy3);
 995             dz33             = _mm_sub_pd(iz3,jz3);
 996
 997             /* Calculate squared distance and things based on it */
 998             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
 999             rsq11            = gmx_mm_calc_rsq_pd(dx11,dy11,dz11);
1000             rsq12            = gmx_mm_calc_rsq_pd(dx12,dy12,dz12);
1001             rsq13            = gmx_mm_calc_rsq_pd(dx13,dy13,dz13);
1002             rsq21            = gmx_mm_calc_rsq_pd(dx21,dy21,dz21);
1003             rsq22            = gmx_mm_calc_rsq_pd(dx22,dy22,dz22);
1004             rsq23            = gmx_mm_calc_rsq_pd(dx23,dy23,dz23);
1005             rsq31            = gmx_mm_calc_rsq_pd(dx31,dy31,dz31);
1006             rsq32            = gmx_mm_calc_rsq_pd(dx32,dy32,dz32);
1007             rsq33            = gmx_mm_calc_rsq_pd(dx33,dy33,dz33);
1008
1009             rinv00           = gmx_mm_invsqrt_pd(rsq00);
1010             rinv11           = gmx_mm_invsqrt_pd(rsq11);
1011             rinv12           = gmx_mm_invsqrt_pd(rsq12);
1012             rinv13           = gmx_mm_invsqrt_pd(rsq13);
1013             rinv21           = gmx_mm_invsqrt_pd(rsq21);
1014             rinv22           = gmx_mm_invsqrt_pd(rsq22);
1015             rinv23           = gmx_mm_invsqrt_pd(rsq23);
1016             rinv31           = gmx_mm_invsqrt_pd(rsq31);
1017             rinv32           = gmx_mm_invsqrt_pd(rsq32);
1018             rinv33           = gmx_mm_invsqrt_pd(rsq33);
1019
1020             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
1021             rinvsq11         = _mm_mul_pd(rinv11,rinv11);
1022             rinvsq12         = _mm_mul_pd(rinv12,rinv12);
1023             rinvsq13         = _mm_mul_pd(rinv13,rinv13);
1024             rinvsq21         = _mm_mul_pd(rinv21,rinv21);
1025             rinvsq22         = _mm_mul_pd(rinv22,rinv22);
1026             rinvsq23         = _mm_mul_pd(rinv23,rinv23);
1027             rinvsq31         = _mm_mul_pd(rinv31,rinv31);
1028             rinvsq32         = _mm_mul_pd(rinv32,rinv32);
1029             rinvsq33         = _mm_mul_pd(rinv33,rinv33);
1030
1031             fjx0             = _mm_setzero_pd();
1032             fjy0             = _mm_setzero_pd();
1033             fjz0             = _mm_setzero_pd();
1034             fjx1             = _mm_setzero_pd();
1035             fjy1             = _mm_setzero_pd();
1036             fjz1             = _mm_setzero_pd();
1037             fjx2             = _mm_setzero_pd();
1038             fjy2             = _mm_setzero_pd();
1039             fjz2             = _mm_setzero_pd();
1040             fjx3             = _mm_setzero_pd();
1041             fjy3             = _mm_setzero_pd();
1042             fjz3             = _mm_setzero_pd();
1043
1044             /**************************
1045              * CALCULATE INTERACTIONS *
1046              **************************/
1047
1048             if (gmx_mm_any_lt(rsq00,rcutoff2))
1049             {
1050
1051             r00              = _mm_mul_pd(rsq00,rinv00);
1052
1053             /* LENNARD-JONES DISPERSION/REPULSION */
1054
1055             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
1056             vvdw6            = _mm_mul_pd(c6_00,rinvsix);
1057             vvdw12           = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
1058             vvdw             = _mm_msub_pd( vvdw12,one_twelfth, _mm_mul_pd(vvdw6,one_sixth) );
1059             fvdw             = _mm_mul_pd(_mm_sub_pd(vvdw12,vvdw6),rinvsq00);
1060
1061             d                = _mm_sub_pd(r00,rswitch);
1062             d                = _mm_max_pd(d,_mm_setzero_pd());
1063             d2               = _mm_mul_pd(d,d);
1064             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
1065
1066             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
1067
1068             /* Evaluate switch function */
1069             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1070             fvdw             = _mm_msub_pd( fvdw,sw , _mm_mul_pd(rinv00,_mm_mul_pd(vvdw,dsw)) );
1071             vvdw             = _mm_mul_pd(vvdw,sw);
1072             cutoff_mask      = _mm_cmplt_pd(rsq00,rcutoff2);
1073
1074             /* Update potential sum for this i atom from the interaction with this j atom. */
1075             vvdw             = _mm_and_pd(vvdw,cutoff_mask);
1076             vvdw             = _mm_unpacklo_pd(vvdw,_mm_setzero_pd());
1077             vvdwsum          = _mm_add_pd(vvdwsum,vvdw);
1078
1079             fscal            = fvdw;
1080
1081             fscal            = _mm_and_pd(fscal,cutoff_mask);
1082
1083             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1084
1085             /* Update vectorial force */
1086             fix0             = _mm_macc_pd(dx00,fscal,fix0);
1087             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
1088             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
1089
1090             fjx0             = _mm_macc_pd(dx00,fscal,fjx0);
1091             fjy0             = _mm_macc_pd(dy00,fscal,fjy0);
1092             fjz0             = _mm_macc_pd(dz00,fscal,fjz0);
1093
1094             }
1095
1096             /**************************
1097              * CALCULATE INTERACTIONS *
1098              **************************/
1099
1100             if (gmx_mm_any_lt(rsq11,rcutoff2))
1101             {
1102
1103             r11              = _mm_mul_pd(rsq11,rinv11);
1104
1105             /* EWALD ELECTROSTATICS */
1106
1107             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1108             ewrt             = _mm_mul_pd(r11,ewtabscale);
1109             ewitab           = _mm_cvttpd_epi32(ewrt);
1110 #ifdef __XOP__
1111             eweps            = _mm_frcz_pd(ewrt);
1112 #else
1113             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1114 #endif
1115             twoeweps         = _mm_add_pd(eweps,eweps);
1116             ewitab           = _mm_slli_epi32(ewitab,2);
1117             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
1118             ewtabD           = _mm_setzero_pd();
1119             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1120             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
1121             ewtabFn          = _mm_setzero_pd();
1122             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1123             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
1124             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
1125             velec            = _mm_mul_pd(qq11,_mm_sub_pd(rinv11,velec));
1126             felec            = _mm_mul_pd(_mm_mul_pd(qq11,rinv11),_mm_sub_pd(rinvsq11,felec));
1127
1128             d                = _mm_sub_pd(r11,rswitch);
1129             d                = _mm_max_pd(d,_mm_setzero_pd());
1130             d2               = _mm_mul_pd(d,d);
1131             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
1132
1133             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
1134
1135             /* Evaluate switch function */
1136             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1137             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv11,_mm_mul_pd(velec,dsw)) );
1138             velec            = _mm_mul_pd(velec,sw);
1139             cutoff_mask      = _mm_cmplt_pd(rsq11,rcutoff2);
1140
1141             /* Update potential sum for this i atom from the interaction with this j atom. */
1142             velec            = _mm_and_pd(velec,cutoff_mask);
1143             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
1144             velecsum         = _mm_add_pd(velecsum,velec);
1145
1146             fscal            = felec;
1147
1148             fscal            = _mm_and_pd(fscal,cutoff_mask);
1149
1150             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1151
1152             /* Update vectorial force */
1153             fix1             = _mm_macc_pd(dx11,fscal,fix1);
1154             fiy1             = _mm_macc_pd(dy11,fscal,fiy1);
1155             fiz1             = _mm_macc_pd(dz11,fscal,fiz1);
1156
1157             fjx1             = _mm_macc_pd(dx11,fscal,fjx1);
1158             fjy1             = _mm_macc_pd(dy11,fscal,fjy1);
1159             fjz1             = _mm_macc_pd(dz11,fscal,fjz1);
1160
1161             }
1162
1163             /**************************
1164              * CALCULATE INTERACTIONS *
1165              **************************/
1166
1167             if (gmx_mm_any_lt(rsq12,rcutoff2))
1168             {
1169
1170             r12              = _mm_mul_pd(rsq12,rinv12);
1171
1172             /* EWALD ELECTROSTATICS */
1173
1174             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1175             ewrt             = _mm_mul_pd(r12,ewtabscale);
1176             ewitab           = _mm_cvttpd_epi32(ewrt);
1177 #ifdef __XOP__
1178             eweps            = _mm_frcz_pd(ewrt);
1179 #else
1180             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1181 #endif
1182             twoeweps         = _mm_add_pd(eweps,eweps);
1183             ewitab           = _mm_slli_epi32(ewitab,2);
1184             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
1185             ewtabD           = _mm_setzero_pd();
1186             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1187             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
1188             ewtabFn          = _mm_setzero_pd();
1189             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1190             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
1191             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
1192             velec            = _mm_mul_pd(qq12,_mm_sub_pd(rinv12,velec));
1193             felec            = _mm_mul_pd(_mm_mul_pd(qq12,rinv12),_mm_sub_pd(rinvsq12,felec));
1194
1195             d                = _mm_sub_pd(r12,rswitch);
1196             d                = _mm_max_pd(d,_mm_setzero_pd());
1197             d2               = _mm_mul_pd(d,d);
1198             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
1199
1200             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
1201
1202             /* Evaluate switch function */
1203             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1204             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv12,_mm_mul_pd(velec,dsw)) );
1205             velec            = _mm_mul_pd(velec,sw);
1206             cutoff_mask      = _mm_cmplt_pd(rsq12,rcutoff2);
1207
1208             /* Update potential sum for this i atom from the interaction with this j atom. */
1209             velec            = _mm_and_pd(velec,cutoff_mask);
1210             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
1211             velecsum         = _mm_add_pd(velecsum,velec);
1212
1213             fscal            = felec;
1214
1215             fscal            = _mm_and_pd(fscal,cutoff_mask);
1216
1217             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1218
1219             /* Update vectorial force */
1220             fix1             = _mm_macc_pd(dx12,fscal,fix1);
1221             fiy1             = _mm_macc_pd(dy12,fscal,fiy1);
1222             fiz1             = _mm_macc_pd(dz12,fscal,fiz1);
1223
1224             fjx2             = _mm_macc_pd(dx12,fscal,fjx2);
1225             fjy2             = _mm_macc_pd(dy12,fscal,fjy2);
1226             fjz2             = _mm_macc_pd(dz12,fscal,fjz2);
1227
1228             }
1229
1230             /**************************
1231              * CALCULATE INTERACTIONS *
1232              **************************/
1233
1234             if (gmx_mm_any_lt(rsq13,rcutoff2))
1235             {
1236
1237             r13              = _mm_mul_pd(rsq13,rinv13);
1238
1239             /* EWALD ELECTROSTATICS */
1240
1241             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1242             ewrt             = _mm_mul_pd(r13,ewtabscale);
1243             ewitab           = _mm_cvttpd_epi32(ewrt);
1244 #ifdef __XOP__
1245             eweps            = _mm_frcz_pd(ewrt);
1246 #else
1247             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1248 #endif
1249             twoeweps         = _mm_add_pd(eweps,eweps);
1250             ewitab           = _mm_slli_epi32(ewitab,2);
1251             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
1252             ewtabD           = _mm_setzero_pd();
1253             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1254             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
1255             ewtabFn          = _mm_setzero_pd();
1256             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1257             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
1258             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
1259             velec            = _mm_mul_pd(qq13,_mm_sub_pd(rinv13,velec));
1260             felec            = _mm_mul_pd(_mm_mul_pd(qq13,rinv13),_mm_sub_pd(rinvsq13,felec));
1261
1262             d                = _mm_sub_pd(r13,rswitch);
1263             d                = _mm_max_pd(d,_mm_setzero_pd());
1264             d2               = _mm_mul_pd(d,d);
1265             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
1266
1267             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
1268
1269             /* Evaluate switch function */
1270             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1271             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv13,_mm_mul_pd(velec,dsw)) );
1272             velec            = _mm_mul_pd(velec,sw);
1273             cutoff_mask      = _mm_cmplt_pd(rsq13,rcutoff2);
1274
1275             /* Update potential sum for this i atom from the interaction with this j atom. */
1276             velec            = _mm_and_pd(velec,cutoff_mask);
1277             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
1278             velecsum         = _mm_add_pd(velecsum,velec);
1279
1280             fscal            = felec;
1281
1282             fscal            = _mm_and_pd(fscal,cutoff_mask);
1283
1284             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1285
1286             /* Update vectorial force */
1287             fix1             = _mm_macc_pd(dx13,fscal,fix1);
1288             fiy1             = _mm_macc_pd(dy13,fscal,fiy1);
1289             fiz1             = _mm_macc_pd(dz13,fscal,fiz1);
1290
1291             fjx3             = _mm_macc_pd(dx13,fscal,fjx3);
1292             fjy3             = _mm_macc_pd(dy13,fscal,fjy3);
1293             fjz3             = _mm_macc_pd(dz13,fscal,fjz3);
1294
1295             }
1296
1297             /**************************
1298              * CALCULATE INTERACTIONS *
1299              **************************/
1300
1301             if (gmx_mm_any_lt(rsq21,rcutoff2))
1302             {
1303
1304             r21              = _mm_mul_pd(rsq21,rinv21);
1305
1306             /* EWALD ELECTROSTATICS */
1307
1308             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1309             ewrt             = _mm_mul_pd(r21,ewtabscale);
1310             ewitab           = _mm_cvttpd_epi32(ewrt);
1311 #ifdef __XOP__
1312             eweps            = _mm_frcz_pd(ewrt);
1313 #else
1314             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1315 #endif
1316             twoeweps         = _mm_add_pd(eweps,eweps);
1317             ewitab           = _mm_slli_epi32(ewitab,2);
1318             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
1319             ewtabD           = _mm_setzero_pd();
1320             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1321             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
1322             ewtabFn          = _mm_setzero_pd();
1323             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1324             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
1325             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
1326             velec            = _mm_mul_pd(qq21,_mm_sub_pd(rinv21,velec));
1327             felec            = _mm_mul_pd(_mm_mul_pd(qq21,rinv21),_mm_sub_pd(rinvsq21,felec));
1328
1329             d                = _mm_sub_pd(r21,rswitch);
1330             d                = _mm_max_pd(d,_mm_setzero_pd());
1331             d2               = _mm_mul_pd(d,d);
1332             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
1333
1334             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
1335
1336             /* Evaluate switch function */
1337             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1338             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv21,_mm_mul_pd(velec,dsw)) );
1339             velec            = _mm_mul_pd(velec,sw);
1340             cutoff_mask      = _mm_cmplt_pd(rsq21,rcutoff2);
1341
1342             /* Update potential sum for this i atom from the interaction with this j atom. */
1343             velec            = _mm_and_pd(velec,cutoff_mask);
1344             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
1345             velecsum         = _mm_add_pd(velecsum,velec);
1346
1347             fscal            = felec;
1348
1349             fscal            = _mm_and_pd(fscal,cutoff_mask);
1350
1351             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1352
1353             /* Update vectorial force */
1354             fix2             = _mm_macc_pd(dx21,fscal,fix2);
1355             fiy2             = _mm_macc_pd(dy21,fscal,fiy2);
1356             fiz2             = _mm_macc_pd(dz21,fscal,fiz2);
1357
1358             fjx1             = _mm_macc_pd(dx21,fscal,fjx1);
1359             fjy1             = _mm_macc_pd(dy21,fscal,fjy1);
1360             fjz1             = _mm_macc_pd(dz21,fscal,fjz1);
1361
1362             }
1363
1364             /**************************
1365              * CALCULATE INTERACTIONS *
1366              **************************/
1367
1368             if (gmx_mm_any_lt(rsq22,rcutoff2))
1369             {
1370
1371             r22              = _mm_mul_pd(rsq22,rinv22);
1372
1373             /* EWALD ELECTROSTATICS */
1374
1375             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1376             ewrt             = _mm_mul_pd(r22,ewtabscale);
1377             ewitab           = _mm_cvttpd_epi32(ewrt);
1378 #ifdef __XOP__
1379             eweps            = _mm_frcz_pd(ewrt);
1380 #else
1381             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1382 #endif
1383             twoeweps         = _mm_add_pd(eweps,eweps);
1384             ewitab           = _mm_slli_epi32(ewitab,2);
1385             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
1386             ewtabD           = _mm_setzero_pd();
1387             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1388             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
1389             ewtabFn          = _mm_setzero_pd();
1390             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1391             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
1392             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
1393             velec            = _mm_mul_pd(qq22,_mm_sub_pd(rinv22,velec));
1394             felec            = _mm_mul_pd(_mm_mul_pd(qq22,rinv22),_mm_sub_pd(rinvsq22,felec));
1395
1396             d                = _mm_sub_pd(r22,rswitch);
1397             d                = _mm_max_pd(d,_mm_setzero_pd());
1398             d2               = _mm_mul_pd(d,d);
1399             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
1400
1401             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
1402
1403             /* Evaluate switch function */
1404             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1405             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv22,_mm_mul_pd(velec,dsw)) );
1406             velec            = _mm_mul_pd(velec,sw);
1407             cutoff_mask      = _mm_cmplt_pd(rsq22,rcutoff2);
1408
1409             /* Update potential sum for this i atom from the interaction with this j atom. */
1410             velec            = _mm_and_pd(velec,cutoff_mask);
1411             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
1412             velecsum         = _mm_add_pd(velecsum,velec);
1413
1414             fscal            = felec;
1415
1416             fscal            = _mm_and_pd(fscal,cutoff_mask);
1417
1418             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1419
1420             /* Update vectorial force */
1421             fix2             = _mm_macc_pd(dx22,fscal,fix2);
1422             fiy2             = _mm_macc_pd(dy22,fscal,fiy2);
1423             fiz2             = _mm_macc_pd(dz22,fscal,fiz2);
1424
1425             fjx2             = _mm_macc_pd(dx22,fscal,fjx2);
1426             fjy2             = _mm_macc_pd(dy22,fscal,fjy2);
1427             fjz2             = _mm_macc_pd(dz22,fscal,fjz2);
1428
1429             }
1430
1431             /**************************
1432              * CALCULATE INTERACTIONS *
1433              **************************/
1434
1435             if (gmx_mm_any_lt(rsq23,rcutoff2))
1436             {
1437
1438             r23              = _mm_mul_pd(rsq23,rinv23);
1439
1440             /* EWALD ELECTROSTATICS */
1441
1442             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1443             ewrt             = _mm_mul_pd(r23,ewtabscale);
1444             ewitab           = _mm_cvttpd_epi32(ewrt);
1445 #ifdef __XOP__
1446             eweps            = _mm_frcz_pd(ewrt);
1447 #else
1448             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1449 #endif
1450             twoeweps         = _mm_add_pd(eweps,eweps);
1451             ewitab           = _mm_slli_epi32(ewitab,2);
1452             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
1453             ewtabD           = _mm_setzero_pd();
1454             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1455             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
1456             ewtabFn          = _mm_setzero_pd();
1457             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1458             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
1459             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
1460             velec            = _mm_mul_pd(qq23,_mm_sub_pd(rinv23,velec));
1461             felec            = _mm_mul_pd(_mm_mul_pd(qq23,rinv23),_mm_sub_pd(rinvsq23,felec));
1462
1463             d                = _mm_sub_pd(r23,rswitch);
1464             d                = _mm_max_pd(d,_mm_setzero_pd());
1465             d2               = _mm_mul_pd(d,d);
1466             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
1467
1468             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
1469
1470             /* Evaluate switch function */
1471             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1472             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv23,_mm_mul_pd(velec,dsw)) );
1473             velec            = _mm_mul_pd(velec,sw);
1474             cutoff_mask      = _mm_cmplt_pd(rsq23,rcutoff2);
1475
1476             /* Update potential sum for this i atom from the interaction with this j atom. */
1477             velec            = _mm_and_pd(velec,cutoff_mask);
1478             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
1479             velecsum         = _mm_add_pd(velecsum,velec);
1480
1481             fscal            = felec;
1482
1483             fscal            = _mm_and_pd(fscal,cutoff_mask);
1484
1485             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1486
1487             /* Update vectorial force */
1488             fix2             = _mm_macc_pd(dx23,fscal,fix2);
1489             fiy2             = _mm_macc_pd(dy23,fscal,fiy2);
1490             fiz2             = _mm_macc_pd(dz23,fscal,fiz2);
1491
1492             fjx3             = _mm_macc_pd(dx23,fscal,fjx3);
1493             fjy3             = _mm_macc_pd(dy23,fscal,fjy3);
1494             fjz3             = _mm_macc_pd(dz23,fscal,fjz3);
1495
1496             }
1497
1498             /**************************
1499              * CALCULATE INTERACTIONS *
1500              **************************/
1501
1502             if (gmx_mm_any_lt(rsq31,rcutoff2))
1503             {
1504
1505             r31              = _mm_mul_pd(rsq31,rinv31);
1506
1507             /* EWALD ELECTROSTATICS */
1508
1509             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1510             ewrt             = _mm_mul_pd(r31,ewtabscale);
1511             ewitab           = _mm_cvttpd_epi32(ewrt);
1512 #ifdef __XOP__
1513             eweps            = _mm_frcz_pd(ewrt);
1514 #else
1515             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1516 #endif
1517             twoeweps         = _mm_add_pd(eweps,eweps);
1518             ewitab           = _mm_slli_epi32(ewitab,2);
1519             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
1520             ewtabD           = _mm_setzero_pd();
1521             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1522             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
1523             ewtabFn          = _mm_setzero_pd();
1524             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1525             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
1526             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
1527             velec            = _mm_mul_pd(qq31,_mm_sub_pd(rinv31,velec));
1528             felec            = _mm_mul_pd(_mm_mul_pd(qq31,rinv31),_mm_sub_pd(rinvsq31,felec));
1529
1530             d                = _mm_sub_pd(r31,rswitch);
1531             d                = _mm_max_pd(d,_mm_setzero_pd());
1532             d2               = _mm_mul_pd(d,d);
1533             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
1534
1535             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
1536
1537             /* Evaluate switch function */
1538             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1539             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv31,_mm_mul_pd(velec,dsw)) );
1540             velec            = _mm_mul_pd(velec,sw);
1541             cutoff_mask      = _mm_cmplt_pd(rsq31,rcutoff2);
1542
1543             /* Update potential sum for this i atom from the interaction with this j atom. */
1544             velec            = _mm_and_pd(velec,cutoff_mask);
1545             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
1546             velecsum         = _mm_add_pd(velecsum,velec);
1547
1548             fscal            = felec;
1549
1550             fscal            = _mm_and_pd(fscal,cutoff_mask);
1551
1552             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1553
1554             /* Update vectorial force */
1555             fix3             = _mm_macc_pd(dx31,fscal,fix3);
1556             fiy3             = _mm_macc_pd(dy31,fscal,fiy3);
1557             fiz3             = _mm_macc_pd(dz31,fscal,fiz3);
1558
1559             fjx1             = _mm_macc_pd(dx31,fscal,fjx1);
1560             fjy1             = _mm_macc_pd(dy31,fscal,fjy1);
1561             fjz1             = _mm_macc_pd(dz31,fscal,fjz1);
1562
1563             }
1564
1565             /**************************
1566              * CALCULATE INTERACTIONS *
1567              **************************/
1568
1569             if (gmx_mm_any_lt(rsq32,rcutoff2))
1570             {
1571
1572             r32              = _mm_mul_pd(rsq32,rinv32);
1573
1574             /* EWALD ELECTROSTATICS */
1575
1576             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1577             ewrt             = _mm_mul_pd(r32,ewtabscale);
1578             ewitab           = _mm_cvttpd_epi32(ewrt);
1579 #ifdef __XOP__
1580             eweps            = _mm_frcz_pd(ewrt);
1581 #else
1582             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1583 #endif
1584             twoeweps         = _mm_add_pd(eweps,eweps);
1585             ewitab           = _mm_slli_epi32(ewitab,2);
1586             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
1587             ewtabD           = _mm_setzero_pd();
1588             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1589             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
1590             ewtabFn          = _mm_setzero_pd();
1591             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1592             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
1593             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
1594             velec            = _mm_mul_pd(qq32,_mm_sub_pd(rinv32,velec));
1595             felec            = _mm_mul_pd(_mm_mul_pd(qq32,rinv32),_mm_sub_pd(rinvsq32,felec));
1596
1597             d                = _mm_sub_pd(r32,rswitch);
1598             d                = _mm_max_pd(d,_mm_setzero_pd());
1599             d2               = _mm_mul_pd(d,d);
1600             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
1601
1602             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
1603
1604             /* Evaluate switch function */
1605             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1606             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv32,_mm_mul_pd(velec,dsw)) );
1607             velec            = _mm_mul_pd(velec,sw);
1608             cutoff_mask      = _mm_cmplt_pd(rsq32,rcutoff2);
1609
1610             /* Update potential sum for this i atom from the interaction with this j atom. */
1611             velec            = _mm_and_pd(velec,cutoff_mask);
1612             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
1613             velecsum         = _mm_add_pd(velecsum,velec);
1614
1615             fscal            = felec;
1616
1617             fscal            = _mm_and_pd(fscal,cutoff_mask);
1618
1619             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1620
1621             /* Update vectorial force */
1622             fix3             = _mm_macc_pd(dx32,fscal,fix3);
1623             fiy3             = _mm_macc_pd(dy32,fscal,fiy3);
1624             fiz3             = _mm_macc_pd(dz32,fscal,fiz3);
1625
1626             fjx2             = _mm_macc_pd(dx32,fscal,fjx2);
1627             fjy2             = _mm_macc_pd(dy32,fscal,fjy2);
1628             fjz2             = _mm_macc_pd(dz32,fscal,fjz2);
1629
1630             }
1631
1632             /**************************
1633              * CALCULATE INTERACTIONS *
1634              **************************/
1635
1636             if (gmx_mm_any_lt(rsq33,rcutoff2))
1637             {
1638
1639             r33              = _mm_mul_pd(rsq33,rinv33);
1640
1641             /* EWALD ELECTROSTATICS */
1642
1643             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1644             ewrt             = _mm_mul_pd(r33,ewtabscale);
1645             ewitab           = _mm_cvttpd_epi32(ewrt);
1646 #ifdef __XOP__
1647             eweps            = _mm_frcz_pd(ewrt);
1648 #else
1649             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1650 #endif
1651             twoeweps         = _mm_add_pd(eweps,eweps);
1652             ewitab           = _mm_slli_epi32(ewitab,2);
1653             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
1654             ewtabD           = _mm_setzero_pd();
1655             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1656             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
1657             ewtabFn          = _mm_setzero_pd();
1658             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1659             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
1660             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
1661             velec            = _mm_mul_pd(qq33,_mm_sub_pd(rinv33,velec));
1662             felec            = _mm_mul_pd(_mm_mul_pd(qq33,rinv33),_mm_sub_pd(rinvsq33,felec));
1663
1664             d                = _mm_sub_pd(r33,rswitch);
1665             d                = _mm_max_pd(d,_mm_setzero_pd());
1666             d2               = _mm_mul_pd(d,d);
1667             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
1668
1669             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
1670
1671             /* Evaluate switch function */
1672             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1673             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv33,_mm_mul_pd(velec,dsw)) );
1674             velec            = _mm_mul_pd(velec,sw);
1675             cutoff_mask      = _mm_cmplt_pd(rsq33,rcutoff2);
1676
1677             /* Update potential sum for this i atom from the interaction with this j atom. */
1678             velec            = _mm_and_pd(velec,cutoff_mask);
1679             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
1680             velecsum         = _mm_add_pd(velecsum,velec);
1681
1682             fscal            = felec;
1683
1684             fscal            = _mm_and_pd(fscal,cutoff_mask);
1685
1686             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1687
1688             /* Update vectorial force */
1689             fix3             = _mm_macc_pd(dx33,fscal,fix3);
1690             fiy3             = _mm_macc_pd(dy33,fscal,fiy3);
1691             fiz3             = _mm_macc_pd(dz33,fscal,fiz3);
1692
1693             fjx3             = _mm_macc_pd(dx33,fscal,fjx3);
1694             fjy3             = _mm_macc_pd(dy33,fscal,fjy3);
1695             fjz3             = _mm_macc_pd(dz33,fscal,fjz3);
1696
1697             }
1698
1699             gmx_mm_decrement_4rvec_1ptr_swizzle_pd(f+j_coord_offsetA,fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,fjx2,fjy2,fjz2,fjx3,fjy3,fjz3);
1700
1701             /* Inner loop uses 677 flops */
1702         }
1703
1704         /* End of innermost loop */
1705
1706         gmx_mm_update_iforce_4atom_swizzle_pd(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
1707                                               f+i_coord_offset,fshift+i_shift_offset);
1708
1709         ggid                        = gid[iidx];
1710         /* Update potential energies */
1711         gmx_mm_update_1pot_pd(velecsum,kernel_data->energygrp_elec+ggid);
1712         gmx_mm_update_1pot_pd(vvdwsum,kernel_data->energygrp_vdw+ggid);
1713
1714         /* Increment number of inner iterations */
1715         inneriter                  += j_index_end - j_index_start;
1716
1717         /* Outer loop uses 26 flops */
1718     }
1719
1720     /* Increment number of outer iterations */
1721     outeriter        += nri;
1722
1723     /* Update outer/inner flops */
1724
1725     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4W4_VF,outeriter*26 + inneriter*677);
1726 }
1727 /*
1728  * Gromacs nonbonded kernel:   nb_kernel_ElecEwSw_VdwLJSw_GeomW4W4_F_avx_128_fma_double
1729  * Electrostatics interaction: Ewald
1730  * VdW interaction:            LennardJones
1731  * Geometry:                   Water4-Water4
1732  * Calculate force/pot:        Force
1733  */
1734 void
1735 nb_kernel_ElecEwSw_VdwLJSw_GeomW4W4_F_avx_128_fma_double
1736                     (t_nblist                    * gmx_restrict       nlist,
1737                      rvec                        * gmx_restrict          xx,
1738                      rvec                        * gmx_restrict          ff,
1739                      t_forcerec                  * gmx_restrict          fr,
1740                      t_mdatoms                   * gmx_restrict     mdatoms,
1741                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
1742                      t_nrnb                      * gmx_restrict        nrnb)
1743 {
1744     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
1745      * just 0 for non-waters.
1746      * Suffixes A,B refer to j loop unrolling done with SSE double precision, e.g. for the two different
1747      * jnr indices corresponding to data put in the four positions in the SIMD register.
1748      */
1749     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
1750     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
1751     int              jnrA,jnrB;
1752     int              j_coord_offsetA,j_coord_offsetB;
1753     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
1754     real             rcutoff_scalar;
1755     real             *shiftvec,*fshift,*x,*f;
1756     __m128d          tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
1757     int              vdwioffset0;
1758     __m128d          ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
1759     int              vdwioffset1;
1760     __m128d          ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
1761     int              vdwioffset2;
1762     __m128d          ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
1763     int              vdwioffset3;
1764     __m128d          ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
1765     int              vdwjidx0A,vdwjidx0B;
1766     __m128d          jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
1767     int              vdwjidx1A,vdwjidx1B;
1768     __m128d          jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
1769     int              vdwjidx2A,vdwjidx2B;
1770     __m128d          jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
1771     int              vdwjidx3A,vdwjidx3B;
1772     __m128d          jx3,jy3,jz3,fjx3,fjy3,fjz3,jq3,isaj3;
1773     __m128d          dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
1774     __m128d          dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11;
1775     __m128d          dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12;
1776     __m128d          dx13,dy13,dz13,rsq13,rinv13,rinvsq13,r13,qq13,c6_13,c12_13;
1777     __m128d          dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21;
1778     __m128d          dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22;
1779     __m128d          dx23,dy23,dz23,rsq23,rinv23,rinvsq23,r23,qq23,c6_23,c12_23;
1780     __m128d          dx31,dy31,dz31,rsq31,rinv31,rinvsq31,r31,qq31,c6_31,c12_31;
1781     __m128d          dx32,dy32,dz32,rsq32,rinv32,rinvsq32,r32,qq32,c6_32,c12_32;
1782     __m128d          dx33,dy33,dz33,rsq33,rinv33,rinvsq33,r33,qq33,c6_33,c12_33;
1783     __m128d          velec,felec,velecsum,facel,crf,krf,krf2;
1784     real             *charge;
1785     int              nvdwtype;
1786     __m128d          rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
1787     int              *vdwtype;
1788     real             *vdwparam;
1789     __m128d          one_sixth   = _mm_set1_pd(1.0/6.0);
1790     __m128d          one_twelfth = _mm_set1_pd(1.0/12.0);
1791     __m128i          ewitab;
1792     __m128d          ewtabscale,eweps,twoeweps,sh_ewald,ewrt,ewtabhalfspace,ewtabF,ewtabFn,ewtabD,ewtabV;
1793     real             *ewtab;
1794     __m128d          rswitch,swV3,swV4,swV5,swF2,swF3,swF4,d,d2,sw,dsw;
1795     real             rswitch_scalar,d_scalar;
1796     __m128d          dummy_mask,cutoff_mask;
1797     __m128d          signbit   = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
1798     __m128d          one     = _mm_set1_pd(1.0);
1799     __m128d          two     = _mm_set1_pd(2.0);
1800     x                = xx[0];
1801     f                = ff[0];
1802
1803     nri              = nlist->nri;
1804     iinr             = nlist->iinr;
1805     jindex           = nlist->jindex;
1806     jjnr             = nlist->jjnr;
1807     shiftidx         = nlist->shift;
1808     gid              = nlist->gid;
1809     shiftvec         = fr->shift_vec[0];
1810     fshift           = fr->fshift[0];
1811     facel            = _mm_set1_pd(fr->epsfac);
1812     charge           = mdatoms->chargeA;
1813     nvdwtype         = fr->ntype;
1814     vdwparam         = fr->nbfp;
1815     vdwtype          = mdatoms->typeA;
1816
1817     sh_ewald         = _mm_set1_pd(fr->ic->sh_ewald);
1818     ewtab            = fr->ic->tabq_coul_FDV0;
1819     ewtabscale       = _mm_set1_pd(fr->ic->tabq_scale);
1820     ewtabhalfspace   = _mm_set1_pd(0.5/fr->ic->tabq_scale);
1821
1822     /* Setup water-specific parameters */
1823     inr              = nlist->iinr[0];
1824     iq1              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+1]));
1825     iq2              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+2]));
1826     iq3              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+3]));
1827     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
1828
1829     jq1              = _mm_set1_pd(charge[inr+1]);
1830     jq2              = _mm_set1_pd(charge[inr+2]);
1831     jq3              = _mm_set1_pd(charge[inr+3]);
1832     vdwjidx0A        = 2*vdwtype[inr+0];
1833     c6_00            = _mm_set1_pd(vdwparam[vdwioffset0+vdwjidx0A]);
1834     c12_00           = _mm_set1_pd(vdwparam[vdwioffset0+vdwjidx0A+1]);
1835     qq11             = _mm_mul_pd(iq1,jq1);
1836     qq12             = _mm_mul_pd(iq1,jq2);
1837     qq13             = _mm_mul_pd(iq1,jq3);
1838     qq21             = _mm_mul_pd(iq2,jq1);
1839     qq22             = _mm_mul_pd(iq2,jq2);
1840     qq23             = _mm_mul_pd(iq2,jq3);
1841     qq31             = _mm_mul_pd(iq3,jq1);
1842     qq32             = _mm_mul_pd(iq3,jq2);
1843     qq33             = _mm_mul_pd(iq3,jq3);
1844
1845     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
1846     rcutoff_scalar   = fr->rcoulomb;
1847     rcutoff          = _mm_set1_pd(rcutoff_scalar);
1848     rcutoff2         = _mm_mul_pd(rcutoff,rcutoff);
1849
1850     rswitch_scalar   = fr->rcoulomb_switch;
1851     rswitch          = _mm_set1_pd(rswitch_scalar);
1852     /* Setup switch parameters */
1853     d_scalar         = rcutoff_scalar-rswitch_scalar;
1854     d                = _mm_set1_pd(d_scalar);
1855     swV3             = _mm_set1_pd(-10.0/(d_scalar*d_scalar*d_scalar));
1856     swV4             = _mm_set1_pd( 15.0/(d_scalar*d_scalar*d_scalar*d_scalar));
1857     swV5             = _mm_set1_pd( -6.0/(d_scalar*d_scalar*d_scalar*d_scalar*d_scalar));
1858     swF2             = _mm_set1_pd(-30.0/(d_scalar*d_scalar*d_scalar));
1859     swF3             = _mm_set1_pd( 60.0/(d_scalar*d_scalar*d_scalar*d_scalar));
1860     swF4             = _mm_set1_pd(-30.0/(d_scalar*d_scalar*d_scalar*d_scalar*d_scalar));
1861
1862     /* Avoid stupid compiler warnings */
1863     jnrA = jnrB = 0;
1864     j_coord_offsetA = 0;
1865     j_coord_offsetB = 0;
1866
1867     outeriter        = 0;
1868     inneriter        = 0;
1869
1870     /* Start outer loop over neighborlists */
1871     for(iidx=0; iidx<nri; iidx++)
1872     {
1873         /* Load shift vector for this list */
1874         i_shift_offset   = DIM*shiftidx[iidx];
1875
1876         /* Load limits for loop over neighbors */
1877         j_index_start    = jindex[iidx];
1878         j_index_end      = jindex[iidx+1];
1879
1880         /* Get outer coordinate index */
1881         inr              = iinr[iidx];
1882         i_coord_offset   = DIM*inr;
1883
1884         /* Load i particle coords and add shift vector */
1885         gmx_mm_load_shift_and_4rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,
1886                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
1887
1888         fix0             = _mm_setzero_pd();
1889         fiy0             = _mm_setzero_pd();
1890         fiz0             = _mm_setzero_pd();
1891         fix1             = _mm_setzero_pd();
1892         fiy1             = _mm_setzero_pd();
1893         fiz1             = _mm_setzero_pd();
1894         fix2             = _mm_setzero_pd();
1895         fiy2             = _mm_setzero_pd();
1896         fiz2             = _mm_setzero_pd();
1897         fix3             = _mm_setzero_pd();
1898         fiy3             = _mm_setzero_pd();
1899         fiz3             = _mm_setzero_pd();
1900
1901         /* Start inner kernel loop */
1902         for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
1903         {
1904
1905             /* Get j neighbor index, and coordinate index */
1906             jnrA             = jjnr[jidx];
1907             jnrB             = jjnr[jidx+1];
1908             j_coord_offsetA  = DIM*jnrA;
1909             j_coord_offsetB  = DIM*jnrB;
1910
1911             /* load j atom coordinates */
1912             gmx_mm_load_4rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
1913                                               &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,
1914                                               &jy2,&jz2,&jx3,&jy3,&jz3);
1915
1916             /* Calculate displacement vector */
1917             dx00             = _mm_sub_pd(ix0,jx0);
1918             dy00             = _mm_sub_pd(iy0,jy0);
1919             dz00             = _mm_sub_pd(iz0,jz0);
1920             dx11             = _mm_sub_pd(ix1,jx1);
1921             dy11             = _mm_sub_pd(iy1,jy1);
1922             dz11             = _mm_sub_pd(iz1,jz1);
1923             dx12             = _mm_sub_pd(ix1,jx2);
1924             dy12             = _mm_sub_pd(iy1,jy2);
1925             dz12             = _mm_sub_pd(iz1,jz2);
1926             dx13             = _mm_sub_pd(ix1,jx3);
1927             dy13             = _mm_sub_pd(iy1,jy3);
1928             dz13             = _mm_sub_pd(iz1,jz3);
1929             dx21             = _mm_sub_pd(ix2,jx1);
1930             dy21             = _mm_sub_pd(iy2,jy1);
1931             dz21             = _mm_sub_pd(iz2,jz1);
1932             dx22             = _mm_sub_pd(ix2,jx2);
1933             dy22             = _mm_sub_pd(iy2,jy2);
1934             dz22             = _mm_sub_pd(iz2,jz2);
1935             dx23             = _mm_sub_pd(ix2,jx3);
1936             dy23             = _mm_sub_pd(iy2,jy3);
1937             dz23             = _mm_sub_pd(iz2,jz3);
1938             dx31             = _mm_sub_pd(ix3,jx1);
1939             dy31             = _mm_sub_pd(iy3,jy1);
1940             dz31             = _mm_sub_pd(iz3,jz1);
1941             dx32             = _mm_sub_pd(ix3,jx2);
1942             dy32             = _mm_sub_pd(iy3,jy2);
1943             dz32             = _mm_sub_pd(iz3,jz2);
1944             dx33             = _mm_sub_pd(ix3,jx3);
1945             dy33             = _mm_sub_pd(iy3,jy3);
1946             dz33             = _mm_sub_pd(iz3,jz3);
1947
1948             /* Calculate squared distance and things based on it */
1949             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
1950             rsq11            = gmx_mm_calc_rsq_pd(dx11,dy11,dz11);
1951             rsq12            = gmx_mm_calc_rsq_pd(dx12,dy12,dz12);
1952             rsq13            = gmx_mm_calc_rsq_pd(dx13,dy13,dz13);
1953             rsq21            = gmx_mm_calc_rsq_pd(dx21,dy21,dz21);
1954             rsq22            = gmx_mm_calc_rsq_pd(dx22,dy22,dz22);
1955             rsq23            = gmx_mm_calc_rsq_pd(dx23,dy23,dz23);
1956             rsq31            = gmx_mm_calc_rsq_pd(dx31,dy31,dz31);
1957             rsq32            = gmx_mm_calc_rsq_pd(dx32,dy32,dz32);
1958             rsq33            = gmx_mm_calc_rsq_pd(dx33,dy33,dz33);
1959
1960             rinv00           = gmx_mm_invsqrt_pd(rsq00);
1961             rinv11           = gmx_mm_invsqrt_pd(rsq11);
1962             rinv12           = gmx_mm_invsqrt_pd(rsq12);
1963             rinv13           = gmx_mm_invsqrt_pd(rsq13);
1964             rinv21           = gmx_mm_invsqrt_pd(rsq21);
1965             rinv22           = gmx_mm_invsqrt_pd(rsq22);
1966             rinv23           = gmx_mm_invsqrt_pd(rsq23);
1967             rinv31           = gmx_mm_invsqrt_pd(rsq31);
1968             rinv32           = gmx_mm_invsqrt_pd(rsq32);
1969             rinv33           = gmx_mm_invsqrt_pd(rsq33);
1970
1971             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
1972             rinvsq11         = _mm_mul_pd(rinv11,rinv11);
1973             rinvsq12         = _mm_mul_pd(rinv12,rinv12);
1974             rinvsq13         = _mm_mul_pd(rinv13,rinv13);
1975             rinvsq21         = _mm_mul_pd(rinv21,rinv21);
1976             rinvsq22         = _mm_mul_pd(rinv22,rinv22);
1977             rinvsq23         = _mm_mul_pd(rinv23,rinv23);
1978             rinvsq31         = _mm_mul_pd(rinv31,rinv31);
1979             rinvsq32         = _mm_mul_pd(rinv32,rinv32);
1980             rinvsq33         = _mm_mul_pd(rinv33,rinv33);
1981
1982             fjx0             = _mm_setzero_pd();
1983             fjy0             = _mm_setzero_pd();
1984             fjz0             = _mm_setzero_pd();
1985             fjx1             = _mm_setzero_pd();
1986             fjy1             = _mm_setzero_pd();
1987             fjz1             = _mm_setzero_pd();
1988             fjx2             = _mm_setzero_pd();
1989             fjy2             = _mm_setzero_pd();
1990             fjz2             = _mm_setzero_pd();
1991             fjx3             = _mm_setzero_pd();
1992             fjy3             = _mm_setzero_pd();
1993             fjz3             = _mm_setzero_pd();
1994
1995             /**************************
1996              * CALCULATE INTERACTIONS *
1997              **************************/
1998
1999             if (gmx_mm_any_lt(rsq00,rcutoff2))
2000             {
2001
2002             r00              = _mm_mul_pd(rsq00,rinv00);
2003
2004             /* LENNARD-JONES DISPERSION/REPULSION */
2005
2006             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
2007             vvdw6            = _mm_mul_pd(c6_00,rinvsix);
2008             vvdw12           = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
2009             vvdw             = _mm_msub_pd( vvdw12,one_twelfth, _mm_mul_pd(vvdw6,one_sixth) );
2010             fvdw             = _mm_mul_pd(_mm_sub_pd(vvdw12,vvdw6),rinvsq00);
2011
2012             d                = _mm_sub_pd(r00,rswitch);
2013             d                = _mm_max_pd(d,_mm_setzero_pd());
2014             d2               = _mm_mul_pd(d,d);
2015             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
2016
2017             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
2018
2019             /* Evaluate switch function */
2020             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2021             fvdw             = _mm_msub_pd( fvdw,sw , _mm_mul_pd(rinv00,_mm_mul_pd(vvdw,dsw)) );
2022             cutoff_mask      = _mm_cmplt_pd(rsq00,rcutoff2);
2023
2024             fscal            = fvdw;
2025
2026             fscal            = _mm_and_pd(fscal,cutoff_mask);
2027
2028             /* Update vectorial force */
2029             fix0             = _mm_macc_pd(dx00,fscal,fix0);
2030             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
2031             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
2032
2033             fjx0             = _mm_macc_pd(dx00,fscal,fjx0);
2034             fjy0             = _mm_macc_pd(dy00,fscal,fjy0);
2035             fjz0             = _mm_macc_pd(dz00,fscal,fjz0);
2036
2037             }
2038
2039             /**************************
2040              * CALCULATE INTERACTIONS *
2041              **************************/
2042
2043             if (gmx_mm_any_lt(rsq11,rcutoff2))
2044             {
2045
2046             r11              = _mm_mul_pd(rsq11,rinv11);
2047
2048             /* EWALD ELECTROSTATICS */
2049
2050             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2051             ewrt             = _mm_mul_pd(r11,ewtabscale);
2052             ewitab           = _mm_cvttpd_epi32(ewrt);
2053 #ifdef __XOP__
2054             eweps            = _mm_frcz_pd(ewrt);
2055 #else
2056             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
2057 #endif
2058             twoeweps         = _mm_add_pd(eweps,eweps);
2059             ewitab           = _mm_slli_epi32(ewitab,2);
2060             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
2061             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
2062             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2063             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
2064             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
2065             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2066             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
2067             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
2068             velec            = _mm_mul_pd(qq11,_mm_sub_pd(rinv11,velec));
2069             felec            = _mm_mul_pd(_mm_mul_pd(qq11,rinv11),_mm_sub_pd(rinvsq11,felec));
2070
2071             d                = _mm_sub_pd(r11,rswitch);
2072             d                = _mm_max_pd(d,_mm_setzero_pd());
2073             d2               = _mm_mul_pd(d,d);
2074             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
2075
2076             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
2077
2078             /* Evaluate switch function */
2079             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2080             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv11,_mm_mul_pd(velec,dsw)) );
2081             cutoff_mask      = _mm_cmplt_pd(rsq11,rcutoff2);
2082
2083             fscal            = felec;
2084
2085             fscal            = _mm_and_pd(fscal,cutoff_mask);
2086
2087             /* Update vectorial force */
2088             fix1             = _mm_macc_pd(dx11,fscal,fix1);
2089             fiy1             = _mm_macc_pd(dy11,fscal,fiy1);
2090             fiz1             = _mm_macc_pd(dz11,fscal,fiz1);
2091
2092             fjx1             = _mm_macc_pd(dx11,fscal,fjx1);
2093             fjy1             = _mm_macc_pd(dy11,fscal,fjy1);
2094             fjz1             = _mm_macc_pd(dz11,fscal,fjz1);
2095
2096             }
2097
2098             /**************************
2099              * CALCULATE INTERACTIONS *
2100              **************************/
2101
2102             if (gmx_mm_any_lt(rsq12,rcutoff2))
2103             {
2104
2105             r12              = _mm_mul_pd(rsq12,rinv12);
2106
2107             /* EWALD ELECTROSTATICS */
2108
2109             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2110             ewrt             = _mm_mul_pd(r12,ewtabscale);
2111             ewitab           = _mm_cvttpd_epi32(ewrt);
2112 #ifdef __XOP__
2113             eweps            = _mm_frcz_pd(ewrt);
2114 #else
2115             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
2116 #endif
2117             twoeweps         = _mm_add_pd(eweps,eweps);
2118             ewitab           = _mm_slli_epi32(ewitab,2);
2119             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
2120             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
2121             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2122             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
2123             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
2124             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2125             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
2126             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
2127             velec            = _mm_mul_pd(qq12,_mm_sub_pd(rinv12,velec));
2128             felec            = _mm_mul_pd(_mm_mul_pd(qq12,rinv12),_mm_sub_pd(rinvsq12,felec));
2129
2130             d                = _mm_sub_pd(r12,rswitch);
2131             d                = _mm_max_pd(d,_mm_setzero_pd());
2132             d2               = _mm_mul_pd(d,d);
2133             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
2134
2135             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
2136
2137             /* Evaluate switch function */
2138             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2139             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv12,_mm_mul_pd(velec,dsw)) );
2140             cutoff_mask      = _mm_cmplt_pd(rsq12,rcutoff2);
2141
2142             fscal            = felec;
2143
2144             fscal            = _mm_and_pd(fscal,cutoff_mask);
2145
2146             /* Update vectorial force */
2147             fix1             = _mm_macc_pd(dx12,fscal,fix1);
2148             fiy1             = _mm_macc_pd(dy12,fscal,fiy1);
2149             fiz1             = _mm_macc_pd(dz12,fscal,fiz1);
2150
2151             fjx2             = _mm_macc_pd(dx12,fscal,fjx2);
2152             fjy2             = _mm_macc_pd(dy12,fscal,fjy2);
2153             fjz2             = _mm_macc_pd(dz12,fscal,fjz2);
2154
2155             }
2156
2157             /**************************
2158              * CALCULATE INTERACTIONS *
2159              **************************/
2160
2161             if (gmx_mm_any_lt(rsq13,rcutoff2))
2162             {
2163
2164             r13              = _mm_mul_pd(rsq13,rinv13);
2165
2166             /* EWALD ELECTROSTATICS */
2167
2168             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2169             ewrt             = _mm_mul_pd(r13,ewtabscale);
2170             ewitab           = _mm_cvttpd_epi32(ewrt);
2171 #ifdef __XOP__
2172             eweps            = _mm_frcz_pd(ewrt);
2173 #else
2174             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
2175 #endif
2176             twoeweps         = _mm_add_pd(eweps,eweps);
2177             ewitab           = _mm_slli_epi32(ewitab,2);
2178             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
2179             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
2180             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2181             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
2182             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
2183             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2184             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
2185             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
2186             velec            = _mm_mul_pd(qq13,_mm_sub_pd(rinv13,velec));
2187             felec            = _mm_mul_pd(_mm_mul_pd(qq13,rinv13),_mm_sub_pd(rinvsq13,felec));
2188
2189             d                = _mm_sub_pd(r13,rswitch);
2190             d                = _mm_max_pd(d,_mm_setzero_pd());
2191             d2               = _mm_mul_pd(d,d);
2192             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
2193
2194             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
2195
2196             /* Evaluate switch function */
2197             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2198             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv13,_mm_mul_pd(velec,dsw)) );
2199             cutoff_mask      = _mm_cmplt_pd(rsq13,rcutoff2);
2200
2201             fscal            = felec;
2202
2203             fscal            = _mm_and_pd(fscal,cutoff_mask);
2204
2205             /* Update vectorial force */
2206             fix1             = _mm_macc_pd(dx13,fscal,fix1);
2207             fiy1             = _mm_macc_pd(dy13,fscal,fiy1);
2208             fiz1             = _mm_macc_pd(dz13,fscal,fiz1);
2209
2210             fjx3             = _mm_macc_pd(dx13,fscal,fjx3);
2211             fjy3             = _mm_macc_pd(dy13,fscal,fjy3);
2212             fjz3             = _mm_macc_pd(dz13,fscal,fjz3);
2213
2214             }
2215
2216             /**************************
2217              * CALCULATE INTERACTIONS *
2218              **************************/
2219
2220             if (gmx_mm_any_lt(rsq21,rcutoff2))
2221             {
2222
2223             r21              = _mm_mul_pd(rsq21,rinv21);
2224
2225             /* EWALD ELECTROSTATICS */
2226
2227             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2228             ewrt             = _mm_mul_pd(r21,ewtabscale);
2229             ewitab           = _mm_cvttpd_epi32(ewrt);
2230 #ifdef __XOP__
2231             eweps            = _mm_frcz_pd(ewrt);
2232 #else
2233             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
2234 #endif
2235             twoeweps         = _mm_add_pd(eweps,eweps);
2236             ewitab           = _mm_slli_epi32(ewitab,2);
2237             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
2238             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
2239             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2240             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
2241             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
2242             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2243             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
2244             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
2245             velec            = _mm_mul_pd(qq21,_mm_sub_pd(rinv21,velec));
2246             felec            = _mm_mul_pd(_mm_mul_pd(qq21,rinv21),_mm_sub_pd(rinvsq21,felec));
2247
2248             d                = _mm_sub_pd(r21,rswitch);
2249             d                = _mm_max_pd(d,_mm_setzero_pd());
2250             d2               = _mm_mul_pd(d,d);
2251             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
2252
2253             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
2254
2255             /* Evaluate switch function */
2256             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2257             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv21,_mm_mul_pd(velec,dsw)) );
2258             cutoff_mask      = _mm_cmplt_pd(rsq21,rcutoff2);
2259
2260             fscal            = felec;
2261
2262             fscal            = _mm_and_pd(fscal,cutoff_mask);
2263
2264             /* Update vectorial force */
2265             fix2             = _mm_macc_pd(dx21,fscal,fix2);
2266             fiy2             = _mm_macc_pd(dy21,fscal,fiy2);
2267             fiz2             = _mm_macc_pd(dz21,fscal,fiz2);
2268
2269             fjx1             = _mm_macc_pd(dx21,fscal,fjx1);
2270             fjy1             = _mm_macc_pd(dy21,fscal,fjy1);
2271             fjz1             = _mm_macc_pd(dz21,fscal,fjz1);
2272
2273             }
2274
2275             /**************************
2276              * CALCULATE INTERACTIONS *
2277              **************************/
2278
2279             if (gmx_mm_any_lt(rsq22,rcutoff2))
2280             {
2281
2282             r22              = _mm_mul_pd(rsq22,rinv22);
2283
2284             /* EWALD ELECTROSTATICS */
2285
2286             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2287             ewrt             = _mm_mul_pd(r22,ewtabscale);
2288             ewitab           = _mm_cvttpd_epi32(ewrt);
2289 #ifdef __XOP__
2290             eweps            = _mm_frcz_pd(ewrt);
2291 #else
2292             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
2293 #endif
2294             twoeweps         = _mm_add_pd(eweps,eweps);
2295             ewitab           = _mm_slli_epi32(ewitab,2);
2296             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
2297             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
2298             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2299             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
2300             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
2301             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2302             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
2303             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
2304             velec            = _mm_mul_pd(qq22,_mm_sub_pd(rinv22,velec));
2305             felec            = _mm_mul_pd(_mm_mul_pd(qq22,rinv22),_mm_sub_pd(rinvsq22,felec));
2306
2307             d                = _mm_sub_pd(r22,rswitch);
2308             d                = _mm_max_pd(d,_mm_setzero_pd());
2309             d2               = _mm_mul_pd(d,d);
2310             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
2311
2312             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
2313
2314             /* Evaluate switch function */
2315             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2316             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv22,_mm_mul_pd(velec,dsw)) );
2317             cutoff_mask      = _mm_cmplt_pd(rsq22,rcutoff2);
2318
2319             fscal            = felec;
2320
2321             fscal            = _mm_and_pd(fscal,cutoff_mask);
2322
2323             /* Update vectorial force */
2324             fix2             = _mm_macc_pd(dx22,fscal,fix2);
2325             fiy2             = _mm_macc_pd(dy22,fscal,fiy2);
2326             fiz2             = _mm_macc_pd(dz22,fscal,fiz2);
2327
2328             fjx2             = _mm_macc_pd(dx22,fscal,fjx2);
2329             fjy2             = _mm_macc_pd(dy22,fscal,fjy2);
2330             fjz2             = _mm_macc_pd(dz22,fscal,fjz2);
2331
2332             }
2333
2334             /**************************
2335              * CALCULATE INTERACTIONS *
2336              **************************/
2337
2338             if (gmx_mm_any_lt(rsq23,rcutoff2))
2339             {
2340
2341             r23              = _mm_mul_pd(rsq23,rinv23);
2342
2343             /* EWALD ELECTROSTATICS */
2344
2345             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2346             ewrt             = _mm_mul_pd(r23,ewtabscale);
2347             ewitab           = _mm_cvttpd_epi32(ewrt);
2348 #ifdef __XOP__
2349             eweps            = _mm_frcz_pd(ewrt);
2350 #else
2351             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
2352 #endif
2353             twoeweps         = _mm_add_pd(eweps,eweps);
2354             ewitab           = _mm_slli_epi32(ewitab,2);
2355             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
2356             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
2357             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2358             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
2359             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
2360             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2361             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
2362             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
2363             velec            = _mm_mul_pd(qq23,_mm_sub_pd(rinv23,velec));
2364             felec            = _mm_mul_pd(_mm_mul_pd(qq23,rinv23),_mm_sub_pd(rinvsq23,felec));
2365
2366             d                = _mm_sub_pd(r23,rswitch);
2367             d                = _mm_max_pd(d,_mm_setzero_pd());
2368             d2               = _mm_mul_pd(d,d);
2369             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
2370
2371             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
2372
2373             /* Evaluate switch function */
2374             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2375             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv23,_mm_mul_pd(velec,dsw)) );
2376             cutoff_mask      = _mm_cmplt_pd(rsq23,rcutoff2);
2377
2378             fscal            = felec;
2379
2380             fscal            = _mm_and_pd(fscal,cutoff_mask);
2381
2382             /* Update vectorial force */
2383             fix2             = _mm_macc_pd(dx23,fscal,fix2);
2384             fiy2             = _mm_macc_pd(dy23,fscal,fiy2);
2385             fiz2             = _mm_macc_pd(dz23,fscal,fiz2);
2386
2387             fjx3             = _mm_macc_pd(dx23,fscal,fjx3);
2388             fjy3             = _mm_macc_pd(dy23,fscal,fjy3);
2389             fjz3             = _mm_macc_pd(dz23,fscal,fjz3);
2390
2391             }
2392
2393             /**************************
2394              * CALCULATE INTERACTIONS *
2395              **************************/
2396
2397             if (gmx_mm_any_lt(rsq31,rcutoff2))
2398             {
2399
2400             r31              = _mm_mul_pd(rsq31,rinv31);
2401
2402             /* EWALD ELECTROSTATICS */
2403
2404             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2405             ewrt             = _mm_mul_pd(r31,ewtabscale);
2406             ewitab           = _mm_cvttpd_epi32(ewrt);
2407 #ifdef __XOP__
2408             eweps            = _mm_frcz_pd(ewrt);
2409 #else
2410             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
2411 #endif
2412             twoeweps         = _mm_add_pd(eweps,eweps);
2413             ewitab           = _mm_slli_epi32(ewitab,2);
2414             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
2415             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
2416             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2417             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
2418             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
2419             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2420             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
2421             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
2422             velec            = _mm_mul_pd(qq31,_mm_sub_pd(rinv31,velec));
2423             felec            = _mm_mul_pd(_mm_mul_pd(qq31,rinv31),_mm_sub_pd(rinvsq31,felec));
2424
2425             d                = _mm_sub_pd(r31,rswitch);
2426             d                = _mm_max_pd(d,_mm_setzero_pd());
2427             d2               = _mm_mul_pd(d,d);
2428             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
2429
2430             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
2431
2432             /* Evaluate switch function */
2433             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2434             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv31,_mm_mul_pd(velec,dsw)) );
2435             cutoff_mask      = _mm_cmplt_pd(rsq31,rcutoff2);
2436
2437             fscal            = felec;
2438
2439             fscal            = _mm_and_pd(fscal,cutoff_mask);
2440
2441             /* Update vectorial force */
2442             fix3             = _mm_macc_pd(dx31,fscal,fix3);
2443             fiy3             = _mm_macc_pd(dy31,fscal,fiy3);
2444             fiz3             = _mm_macc_pd(dz31,fscal,fiz3);
2445
2446             fjx1             = _mm_macc_pd(dx31,fscal,fjx1);
2447             fjy1             = _mm_macc_pd(dy31,fscal,fjy1);
2448             fjz1             = _mm_macc_pd(dz31,fscal,fjz1);
2449
2450             }
2451
2452             /**************************
2453              * CALCULATE INTERACTIONS *
2454              **************************/
2455
2456             if (gmx_mm_any_lt(rsq32,rcutoff2))
2457             {
2458
2459             r32              = _mm_mul_pd(rsq32,rinv32);
2460
2461             /* EWALD ELECTROSTATICS */
2462
2463             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2464             ewrt             = _mm_mul_pd(r32,ewtabscale);
2465             ewitab           = _mm_cvttpd_epi32(ewrt);
2466 #ifdef __XOP__
2467             eweps            = _mm_frcz_pd(ewrt);
2468 #else
2469             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
2470 #endif
2471             twoeweps         = _mm_add_pd(eweps,eweps);
2472             ewitab           = _mm_slli_epi32(ewitab,2);
2473             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
2474             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
2475             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2476             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
2477             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
2478             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2479             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
2480             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
2481             velec            = _mm_mul_pd(qq32,_mm_sub_pd(rinv32,velec));
2482             felec            = _mm_mul_pd(_mm_mul_pd(qq32,rinv32),_mm_sub_pd(rinvsq32,felec));
2483
2484             d                = _mm_sub_pd(r32,rswitch);
2485             d                = _mm_max_pd(d,_mm_setzero_pd());
2486             d2               = _mm_mul_pd(d,d);
2487             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
2488
2489             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
2490
2491             /* Evaluate switch function */
2492             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2493             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv32,_mm_mul_pd(velec,dsw)) );
2494             cutoff_mask      = _mm_cmplt_pd(rsq32,rcutoff2);
2495
2496             fscal            = felec;
2497
2498             fscal            = _mm_and_pd(fscal,cutoff_mask);
2499
2500             /* Update vectorial force */
2501             fix3             = _mm_macc_pd(dx32,fscal,fix3);
2502             fiy3             = _mm_macc_pd(dy32,fscal,fiy3);
2503             fiz3             = _mm_macc_pd(dz32,fscal,fiz3);
2504
2505             fjx2             = _mm_macc_pd(dx32,fscal,fjx2);
2506             fjy2             = _mm_macc_pd(dy32,fscal,fjy2);
2507             fjz2             = _mm_macc_pd(dz32,fscal,fjz2);
2508
2509             }
2510
2511             /**************************
2512              * CALCULATE INTERACTIONS *
2513              **************************/
2514
2515             if (gmx_mm_any_lt(rsq33,rcutoff2))
2516             {
2517
2518             r33              = _mm_mul_pd(rsq33,rinv33);
2519
2520             /* EWALD ELECTROSTATICS */
2521
2522             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2523             ewrt             = _mm_mul_pd(r33,ewtabscale);
2524             ewitab           = _mm_cvttpd_epi32(ewrt);
2525 #ifdef __XOP__
2526             eweps            = _mm_frcz_pd(ewrt);
2527 #else
2528             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
2529 #endif
2530             twoeweps         = _mm_add_pd(eweps,eweps);
2531             ewitab           = _mm_slli_epi32(ewitab,2);
2532             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
2533             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
2534             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2535             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
2536             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
2537             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2538             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
2539             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
2540             velec            = _mm_mul_pd(qq33,_mm_sub_pd(rinv33,velec));
2541             felec            = _mm_mul_pd(_mm_mul_pd(qq33,rinv33),_mm_sub_pd(rinvsq33,felec));
2542
2543             d                = _mm_sub_pd(r33,rswitch);
2544             d                = _mm_max_pd(d,_mm_setzero_pd());
2545             d2               = _mm_mul_pd(d,d);
2546             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
2547
2548             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
2549
2550             /* Evaluate switch function */
2551             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2552             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv33,_mm_mul_pd(velec,dsw)) );
2553             cutoff_mask      = _mm_cmplt_pd(rsq33,rcutoff2);
2554
2555             fscal            = felec;
2556
2557             fscal            = _mm_and_pd(fscal,cutoff_mask);
2558
2559             /* Update vectorial force */
2560             fix3             = _mm_macc_pd(dx33,fscal,fix3);
2561             fiy3             = _mm_macc_pd(dy33,fscal,fiy3);
2562             fiz3             = _mm_macc_pd(dz33,fscal,fiz3);
2563
2564             fjx3             = _mm_macc_pd(dx33,fscal,fjx3);
2565             fjy3             = _mm_macc_pd(dy33,fscal,fjy3);
2566             fjz3             = _mm_macc_pd(dz33,fscal,fjz3);
2567
2568             }
2569
2570             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);
2571
2572             /* Inner loop uses 647 flops */
2573         }
2574
2575         if(jidx<j_index_end)
2576         {
2577
2578             jnrA             = jjnr[jidx];
2579             j_coord_offsetA  = DIM*jnrA;
2580
2581             /* load j atom coordinates */
2582             gmx_mm_load_4rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
2583                                               &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,
2584                                               &jy2,&jz2,&jx3,&jy3,&jz3);
2585
2586             /* Calculate displacement vector */
2587             dx00             = _mm_sub_pd(ix0,jx0);
2588             dy00             = _mm_sub_pd(iy0,jy0);
2589             dz00             = _mm_sub_pd(iz0,jz0);
2590             dx11             = _mm_sub_pd(ix1,jx1);
2591             dy11             = _mm_sub_pd(iy1,jy1);
2592             dz11             = _mm_sub_pd(iz1,jz1);
2593             dx12             = _mm_sub_pd(ix1,jx2);
2594             dy12             = _mm_sub_pd(iy1,jy2);
2595             dz12             = _mm_sub_pd(iz1,jz2);
2596             dx13             = _mm_sub_pd(ix1,jx3);
2597             dy13             = _mm_sub_pd(iy1,jy3);
2598             dz13             = _mm_sub_pd(iz1,jz3);
2599             dx21             = _mm_sub_pd(ix2,jx1);
2600             dy21             = _mm_sub_pd(iy2,jy1);
2601             dz21             = _mm_sub_pd(iz2,jz1);
2602             dx22             = _mm_sub_pd(ix2,jx2);
2603             dy22             = _mm_sub_pd(iy2,jy2);
2604             dz22             = _mm_sub_pd(iz2,jz2);
2605             dx23             = _mm_sub_pd(ix2,jx3);
2606             dy23             = _mm_sub_pd(iy2,jy3);
2607             dz23             = _mm_sub_pd(iz2,jz3);
2608             dx31             = _mm_sub_pd(ix3,jx1);
2609             dy31             = _mm_sub_pd(iy3,jy1);
2610             dz31             = _mm_sub_pd(iz3,jz1);
2611             dx32             = _mm_sub_pd(ix3,jx2);
2612             dy32             = _mm_sub_pd(iy3,jy2);
2613             dz32             = _mm_sub_pd(iz3,jz2);
2614             dx33             = _mm_sub_pd(ix3,jx3);
2615             dy33             = _mm_sub_pd(iy3,jy3);
2616             dz33             = _mm_sub_pd(iz3,jz3);
2617
2618             /* Calculate squared distance and things based on it */
2619             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
2620             rsq11            = gmx_mm_calc_rsq_pd(dx11,dy11,dz11);
2621             rsq12            = gmx_mm_calc_rsq_pd(dx12,dy12,dz12);
2622             rsq13            = gmx_mm_calc_rsq_pd(dx13,dy13,dz13);
2623             rsq21            = gmx_mm_calc_rsq_pd(dx21,dy21,dz21);
2624             rsq22            = gmx_mm_calc_rsq_pd(dx22,dy22,dz22);
2625             rsq23            = gmx_mm_calc_rsq_pd(dx23,dy23,dz23);
2626             rsq31            = gmx_mm_calc_rsq_pd(dx31,dy31,dz31);
2627             rsq32            = gmx_mm_calc_rsq_pd(dx32,dy32,dz32);
2628             rsq33            = gmx_mm_calc_rsq_pd(dx33,dy33,dz33);
2629
2630             rinv00           = gmx_mm_invsqrt_pd(rsq00);
2631             rinv11           = gmx_mm_invsqrt_pd(rsq11);
2632             rinv12           = gmx_mm_invsqrt_pd(rsq12);
2633             rinv13           = gmx_mm_invsqrt_pd(rsq13);
2634             rinv21           = gmx_mm_invsqrt_pd(rsq21);
2635             rinv22           = gmx_mm_invsqrt_pd(rsq22);
2636             rinv23           = gmx_mm_invsqrt_pd(rsq23);
2637             rinv31           = gmx_mm_invsqrt_pd(rsq31);
2638             rinv32           = gmx_mm_invsqrt_pd(rsq32);
2639             rinv33           = gmx_mm_invsqrt_pd(rsq33);
2640
2641             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
2642             rinvsq11         = _mm_mul_pd(rinv11,rinv11);
2643             rinvsq12         = _mm_mul_pd(rinv12,rinv12);
2644             rinvsq13         = _mm_mul_pd(rinv13,rinv13);
2645             rinvsq21         = _mm_mul_pd(rinv21,rinv21);
2646             rinvsq22         = _mm_mul_pd(rinv22,rinv22);
2647             rinvsq23         = _mm_mul_pd(rinv23,rinv23);
2648             rinvsq31         = _mm_mul_pd(rinv31,rinv31);
2649             rinvsq32         = _mm_mul_pd(rinv32,rinv32);
2650             rinvsq33         = _mm_mul_pd(rinv33,rinv33);
2651
2652             fjx0             = _mm_setzero_pd();
2653             fjy0             = _mm_setzero_pd();
2654             fjz0             = _mm_setzero_pd();
2655             fjx1             = _mm_setzero_pd();
2656             fjy1             = _mm_setzero_pd();
2657             fjz1             = _mm_setzero_pd();
2658             fjx2             = _mm_setzero_pd();
2659             fjy2             = _mm_setzero_pd();
2660             fjz2             = _mm_setzero_pd();
2661             fjx3             = _mm_setzero_pd();
2662             fjy3             = _mm_setzero_pd();
2663             fjz3             = _mm_setzero_pd();
2664
2665             /**************************
2666              * CALCULATE INTERACTIONS *
2667              **************************/
2668
2669             if (gmx_mm_any_lt(rsq00,rcutoff2))
2670             {
2671
2672             r00              = _mm_mul_pd(rsq00,rinv00);
2673
2674             /* LENNARD-JONES DISPERSION/REPULSION */
2675
2676             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
2677             vvdw6            = _mm_mul_pd(c6_00,rinvsix);
2678             vvdw12           = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
2679             vvdw             = _mm_msub_pd( vvdw12,one_twelfth, _mm_mul_pd(vvdw6,one_sixth) );
2680             fvdw             = _mm_mul_pd(_mm_sub_pd(vvdw12,vvdw6),rinvsq00);
2681
2682             d                = _mm_sub_pd(r00,rswitch);
2683             d                = _mm_max_pd(d,_mm_setzero_pd());
2684             d2               = _mm_mul_pd(d,d);
2685             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
2686
2687             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
2688
2689             /* Evaluate switch function */
2690             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2691             fvdw             = _mm_msub_pd( fvdw,sw , _mm_mul_pd(rinv00,_mm_mul_pd(vvdw,dsw)) );
2692             cutoff_mask      = _mm_cmplt_pd(rsq00,rcutoff2);
2693
2694             fscal            = fvdw;
2695
2696             fscal            = _mm_and_pd(fscal,cutoff_mask);
2697
2698             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
2699
2700             /* Update vectorial force */
2701             fix0             = _mm_macc_pd(dx00,fscal,fix0);
2702             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
2703             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
2704
2705             fjx0             = _mm_macc_pd(dx00,fscal,fjx0);
2706             fjy0             = _mm_macc_pd(dy00,fscal,fjy0);
2707             fjz0             = _mm_macc_pd(dz00,fscal,fjz0);
2708
2709             }
2710
2711             /**************************
2712              * CALCULATE INTERACTIONS *
2713              **************************/
2714
2715             if (gmx_mm_any_lt(rsq11,rcutoff2))
2716             {
2717
2718             r11              = _mm_mul_pd(rsq11,rinv11);
2719
2720             /* EWALD ELECTROSTATICS */
2721
2722             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2723             ewrt             = _mm_mul_pd(r11,ewtabscale);
2724             ewitab           = _mm_cvttpd_epi32(ewrt);
2725 #ifdef __XOP__
2726             eweps            = _mm_frcz_pd(ewrt);
2727 #else
2728             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
2729 #endif
2730             twoeweps         = _mm_add_pd(eweps,eweps);
2731             ewitab           = _mm_slli_epi32(ewitab,2);
2732             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
2733             ewtabD           = _mm_setzero_pd();
2734             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2735             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
2736             ewtabFn          = _mm_setzero_pd();
2737             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2738             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
2739             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
2740             velec            = _mm_mul_pd(qq11,_mm_sub_pd(rinv11,velec));
2741             felec            = _mm_mul_pd(_mm_mul_pd(qq11,rinv11),_mm_sub_pd(rinvsq11,felec));
2742
2743             d                = _mm_sub_pd(r11,rswitch);
2744             d                = _mm_max_pd(d,_mm_setzero_pd());
2745             d2               = _mm_mul_pd(d,d);
2746             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
2747
2748             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
2749
2750             /* Evaluate switch function */
2751             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2752             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv11,_mm_mul_pd(velec,dsw)) );
2753             cutoff_mask      = _mm_cmplt_pd(rsq11,rcutoff2);
2754
2755             fscal            = felec;
2756
2757             fscal            = _mm_and_pd(fscal,cutoff_mask);
2758
2759             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
2760
2761             /* Update vectorial force */
2762             fix1             = _mm_macc_pd(dx11,fscal,fix1);
2763             fiy1             = _mm_macc_pd(dy11,fscal,fiy1);
2764             fiz1             = _mm_macc_pd(dz11,fscal,fiz1);
2765
2766             fjx1             = _mm_macc_pd(dx11,fscal,fjx1);
2767             fjy1             = _mm_macc_pd(dy11,fscal,fjy1);
2768             fjz1             = _mm_macc_pd(dz11,fscal,fjz1);
2769
2770             }
2771
2772             /**************************
2773              * CALCULATE INTERACTIONS *
2774              **************************/
2775
2776             if (gmx_mm_any_lt(rsq12,rcutoff2))
2777             {
2778
2779             r12              = _mm_mul_pd(rsq12,rinv12);
2780
2781             /* EWALD ELECTROSTATICS */
2782
2783             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2784             ewrt             = _mm_mul_pd(r12,ewtabscale);
2785             ewitab           = _mm_cvttpd_epi32(ewrt);
2786 #ifdef __XOP__
2787             eweps            = _mm_frcz_pd(ewrt);
2788 #else
2789             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
2790 #endif
2791             twoeweps         = _mm_add_pd(eweps,eweps);
2792             ewitab           = _mm_slli_epi32(ewitab,2);
2793             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
2794             ewtabD           = _mm_setzero_pd();
2795             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2796             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
2797             ewtabFn          = _mm_setzero_pd();
2798             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2799             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
2800             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
2801             velec            = _mm_mul_pd(qq12,_mm_sub_pd(rinv12,velec));
2802             felec            = _mm_mul_pd(_mm_mul_pd(qq12,rinv12),_mm_sub_pd(rinvsq12,felec));
2803
2804             d                = _mm_sub_pd(r12,rswitch);
2805             d                = _mm_max_pd(d,_mm_setzero_pd());
2806             d2               = _mm_mul_pd(d,d);
2807             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
2808
2809             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
2810
2811             /* Evaluate switch function */
2812             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2813             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv12,_mm_mul_pd(velec,dsw)) );
2814             cutoff_mask      = _mm_cmplt_pd(rsq12,rcutoff2);
2815
2816             fscal            = felec;
2817
2818             fscal            = _mm_and_pd(fscal,cutoff_mask);
2819
2820             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
2821
2822             /* Update vectorial force */
2823             fix1             = _mm_macc_pd(dx12,fscal,fix1);
2824             fiy1             = _mm_macc_pd(dy12,fscal,fiy1);
2825             fiz1             = _mm_macc_pd(dz12,fscal,fiz1);
2826
2827             fjx2             = _mm_macc_pd(dx12,fscal,fjx2);
2828             fjy2             = _mm_macc_pd(dy12,fscal,fjy2);
2829             fjz2             = _mm_macc_pd(dz12,fscal,fjz2);
2830
2831             }
2832
2833             /**************************
2834              * CALCULATE INTERACTIONS *
2835              **************************/
2836
2837             if (gmx_mm_any_lt(rsq13,rcutoff2))
2838             {
2839
2840             r13              = _mm_mul_pd(rsq13,rinv13);
2841
2842             /* EWALD ELECTROSTATICS */
2843
2844             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2845             ewrt             = _mm_mul_pd(r13,ewtabscale);
2846             ewitab           = _mm_cvttpd_epi32(ewrt);
2847 #ifdef __XOP__
2848             eweps            = _mm_frcz_pd(ewrt);
2849 #else
2850             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
2851 #endif
2852             twoeweps         = _mm_add_pd(eweps,eweps);
2853             ewitab           = _mm_slli_epi32(ewitab,2);
2854             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
2855             ewtabD           = _mm_setzero_pd();
2856             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2857             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
2858             ewtabFn          = _mm_setzero_pd();
2859             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2860             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
2861             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
2862             velec            = _mm_mul_pd(qq13,_mm_sub_pd(rinv13,velec));
2863             felec            = _mm_mul_pd(_mm_mul_pd(qq13,rinv13),_mm_sub_pd(rinvsq13,felec));
2864
2865             d                = _mm_sub_pd(r13,rswitch);
2866             d                = _mm_max_pd(d,_mm_setzero_pd());
2867             d2               = _mm_mul_pd(d,d);
2868             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
2869
2870             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
2871
2872             /* Evaluate switch function */
2873             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2874             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv13,_mm_mul_pd(velec,dsw)) );
2875             cutoff_mask      = _mm_cmplt_pd(rsq13,rcutoff2);
2876
2877             fscal            = felec;
2878
2879             fscal            = _mm_and_pd(fscal,cutoff_mask);
2880
2881             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
2882
2883             /* Update vectorial force */
2884             fix1             = _mm_macc_pd(dx13,fscal,fix1);
2885             fiy1             = _mm_macc_pd(dy13,fscal,fiy1);
2886             fiz1             = _mm_macc_pd(dz13,fscal,fiz1);
2887
2888             fjx3             = _mm_macc_pd(dx13,fscal,fjx3);
2889             fjy3             = _mm_macc_pd(dy13,fscal,fjy3);
2890             fjz3             = _mm_macc_pd(dz13,fscal,fjz3);
2891
2892             }
2893
2894             /**************************
2895              * CALCULATE INTERACTIONS *
2896              **************************/
2897
2898             if (gmx_mm_any_lt(rsq21,rcutoff2))
2899             {
2900
2901             r21              = _mm_mul_pd(rsq21,rinv21);
2902
2903             /* EWALD ELECTROSTATICS */
2904
2905             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2906             ewrt             = _mm_mul_pd(r21,ewtabscale);
2907             ewitab           = _mm_cvttpd_epi32(ewrt);
2908 #ifdef __XOP__
2909             eweps            = _mm_frcz_pd(ewrt);
2910 #else
2911             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
2912 #endif
2913             twoeweps         = _mm_add_pd(eweps,eweps);
2914             ewitab           = _mm_slli_epi32(ewitab,2);
2915             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
2916             ewtabD           = _mm_setzero_pd();
2917             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2918             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
2919             ewtabFn          = _mm_setzero_pd();
2920             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2921             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
2922             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
2923             velec            = _mm_mul_pd(qq21,_mm_sub_pd(rinv21,velec));
2924             felec            = _mm_mul_pd(_mm_mul_pd(qq21,rinv21),_mm_sub_pd(rinvsq21,felec));
2925
2926             d                = _mm_sub_pd(r21,rswitch);
2927             d                = _mm_max_pd(d,_mm_setzero_pd());
2928             d2               = _mm_mul_pd(d,d);
2929             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
2930
2931             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
2932
2933             /* Evaluate switch function */
2934             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2935             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv21,_mm_mul_pd(velec,dsw)) );
2936             cutoff_mask      = _mm_cmplt_pd(rsq21,rcutoff2);
2937
2938             fscal            = felec;
2939
2940             fscal            = _mm_and_pd(fscal,cutoff_mask);
2941
2942             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
2943
2944             /* Update vectorial force */
2945             fix2             = _mm_macc_pd(dx21,fscal,fix2);
2946             fiy2             = _mm_macc_pd(dy21,fscal,fiy2);
2947             fiz2             = _mm_macc_pd(dz21,fscal,fiz2);
2948
2949             fjx1             = _mm_macc_pd(dx21,fscal,fjx1);
2950             fjy1             = _mm_macc_pd(dy21,fscal,fjy1);
2951             fjz1             = _mm_macc_pd(dz21,fscal,fjz1);
2952
2953             }
2954
2955             /**************************
2956              * CALCULATE INTERACTIONS *
2957              **************************/
2958
2959             if (gmx_mm_any_lt(rsq22,rcutoff2))
2960             {
2961
2962             r22              = _mm_mul_pd(rsq22,rinv22);
2963
2964             /* EWALD ELECTROSTATICS */
2965
2966             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2967             ewrt             = _mm_mul_pd(r22,ewtabscale);
2968             ewitab           = _mm_cvttpd_epi32(ewrt);
2969 #ifdef __XOP__
2970             eweps            = _mm_frcz_pd(ewrt);
2971 #else
2972             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
2973 #endif
2974             twoeweps         = _mm_add_pd(eweps,eweps);
2975             ewitab           = _mm_slli_epi32(ewitab,2);
2976             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
2977             ewtabD           = _mm_setzero_pd();
2978             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2979             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
2980             ewtabFn          = _mm_setzero_pd();
2981             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2982             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
2983             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
2984             velec            = _mm_mul_pd(qq22,_mm_sub_pd(rinv22,velec));
2985             felec            = _mm_mul_pd(_mm_mul_pd(qq22,rinv22),_mm_sub_pd(rinvsq22,felec));
2986
2987             d                = _mm_sub_pd(r22,rswitch);
2988             d                = _mm_max_pd(d,_mm_setzero_pd());
2989             d2               = _mm_mul_pd(d,d);
2990             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
2991
2992             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
2993
2994             /* Evaluate switch function */
2995             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2996             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv22,_mm_mul_pd(velec,dsw)) );
2997             cutoff_mask      = _mm_cmplt_pd(rsq22,rcutoff2);
2998
2999             fscal            = felec;
3000
3001             fscal            = _mm_and_pd(fscal,cutoff_mask);
3002
3003             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
3004
3005             /* Update vectorial force */
3006             fix2             = _mm_macc_pd(dx22,fscal,fix2);
3007             fiy2             = _mm_macc_pd(dy22,fscal,fiy2);
3008             fiz2             = _mm_macc_pd(dz22,fscal,fiz2);
3009
3010             fjx2             = _mm_macc_pd(dx22,fscal,fjx2);
3011             fjy2             = _mm_macc_pd(dy22,fscal,fjy2);
3012             fjz2             = _mm_macc_pd(dz22,fscal,fjz2);
3013
3014             }
3015
3016             /**************************
3017              * CALCULATE INTERACTIONS *
3018              **************************/
3019
3020             if (gmx_mm_any_lt(rsq23,rcutoff2))
3021             {
3022
3023             r23              = _mm_mul_pd(rsq23,rinv23);
3024
3025             /* EWALD ELECTROSTATICS */
3026
3027             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
3028             ewrt             = _mm_mul_pd(r23,ewtabscale);
3029             ewitab           = _mm_cvttpd_epi32(ewrt);
3030 #ifdef __XOP__
3031             eweps            = _mm_frcz_pd(ewrt);
3032 #else
3033             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
3034 #endif
3035             twoeweps         = _mm_add_pd(eweps,eweps);
3036             ewitab           = _mm_slli_epi32(ewitab,2);
3037             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
3038             ewtabD           = _mm_setzero_pd();
3039             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
3040             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
3041             ewtabFn          = _mm_setzero_pd();
3042             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
3043             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
3044             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
3045             velec            = _mm_mul_pd(qq23,_mm_sub_pd(rinv23,velec));
3046             felec            = _mm_mul_pd(_mm_mul_pd(qq23,rinv23),_mm_sub_pd(rinvsq23,felec));
3047
3048             d                = _mm_sub_pd(r23,rswitch);
3049             d                = _mm_max_pd(d,_mm_setzero_pd());
3050             d2               = _mm_mul_pd(d,d);
3051             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
3052
3053             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
3054
3055             /* Evaluate switch function */
3056             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
3057             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv23,_mm_mul_pd(velec,dsw)) );
3058             cutoff_mask      = _mm_cmplt_pd(rsq23,rcutoff2);
3059
3060             fscal            = felec;
3061
3062             fscal            = _mm_and_pd(fscal,cutoff_mask);
3063
3064             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
3065
3066             /* Update vectorial force */
3067             fix2             = _mm_macc_pd(dx23,fscal,fix2);
3068             fiy2             = _mm_macc_pd(dy23,fscal,fiy2);
3069             fiz2             = _mm_macc_pd(dz23,fscal,fiz2);
3070
3071             fjx3             = _mm_macc_pd(dx23,fscal,fjx3);
3072             fjy3             = _mm_macc_pd(dy23,fscal,fjy3);
3073             fjz3             = _mm_macc_pd(dz23,fscal,fjz3);
3074
3075             }
3076
3077             /**************************
3078              * CALCULATE INTERACTIONS *
3079              **************************/
3080
3081             if (gmx_mm_any_lt(rsq31,rcutoff2))
3082             {
3083
3084             r31              = _mm_mul_pd(rsq31,rinv31);
3085
3086             /* EWALD ELECTROSTATICS */
3087
3088             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
3089             ewrt             = _mm_mul_pd(r31,ewtabscale);
3090             ewitab           = _mm_cvttpd_epi32(ewrt);
3091 #ifdef __XOP__
3092             eweps            = _mm_frcz_pd(ewrt);
3093 #else
3094             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
3095 #endif
3096             twoeweps         = _mm_add_pd(eweps,eweps);
3097             ewitab           = _mm_slli_epi32(ewitab,2);
3098             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
3099             ewtabD           = _mm_setzero_pd();
3100             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
3101             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
3102             ewtabFn          = _mm_setzero_pd();
3103             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
3104             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
3105             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
3106             velec            = _mm_mul_pd(qq31,_mm_sub_pd(rinv31,velec));
3107             felec            = _mm_mul_pd(_mm_mul_pd(qq31,rinv31),_mm_sub_pd(rinvsq31,felec));
3108
3109             d                = _mm_sub_pd(r31,rswitch);
3110             d                = _mm_max_pd(d,_mm_setzero_pd());
3111             d2               = _mm_mul_pd(d,d);
3112             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
3113
3114             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
3115
3116             /* Evaluate switch function */
3117             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
3118             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv31,_mm_mul_pd(velec,dsw)) );
3119             cutoff_mask      = _mm_cmplt_pd(rsq31,rcutoff2);
3120
3121             fscal            = felec;
3122
3123             fscal            = _mm_and_pd(fscal,cutoff_mask);
3124
3125             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
3126
3127             /* Update vectorial force */
3128             fix3             = _mm_macc_pd(dx31,fscal,fix3);
3129             fiy3             = _mm_macc_pd(dy31,fscal,fiy3);
3130             fiz3             = _mm_macc_pd(dz31,fscal,fiz3);
3131
3132             fjx1             = _mm_macc_pd(dx31,fscal,fjx1);
3133             fjy1             = _mm_macc_pd(dy31,fscal,fjy1);
3134             fjz1             = _mm_macc_pd(dz31,fscal,fjz1);
3135
3136             }
3137
3138             /**************************
3139              * CALCULATE INTERACTIONS *
3140              **************************/
3141
3142             if (gmx_mm_any_lt(rsq32,rcutoff2))
3143             {
3144
3145             r32              = _mm_mul_pd(rsq32,rinv32);
3146
3147             /* EWALD ELECTROSTATICS */
3148
3149             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
3150             ewrt             = _mm_mul_pd(r32,ewtabscale);
3151             ewitab           = _mm_cvttpd_epi32(ewrt);
3152 #ifdef __XOP__
3153             eweps            = _mm_frcz_pd(ewrt);
3154 #else
3155             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
3156 #endif
3157             twoeweps         = _mm_add_pd(eweps,eweps);
3158             ewitab           = _mm_slli_epi32(ewitab,2);
3159             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
3160             ewtabD           = _mm_setzero_pd();
3161             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
3162             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
3163             ewtabFn          = _mm_setzero_pd();
3164             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
3165             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
3166             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
3167             velec            = _mm_mul_pd(qq32,_mm_sub_pd(rinv32,velec));
3168             felec            = _mm_mul_pd(_mm_mul_pd(qq32,rinv32),_mm_sub_pd(rinvsq32,felec));
3169
3170             d                = _mm_sub_pd(r32,rswitch);
3171             d                = _mm_max_pd(d,_mm_setzero_pd());
3172             d2               = _mm_mul_pd(d,d);
3173             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
3174
3175             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
3176
3177             /* Evaluate switch function */
3178             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
3179             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv32,_mm_mul_pd(velec,dsw)) );
3180             cutoff_mask      = _mm_cmplt_pd(rsq32,rcutoff2);
3181
3182             fscal            = felec;
3183
3184             fscal            = _mm_and_pd(fscal,cutoff_mask);
3185
3186             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
3187
3188             /* Update vectorial force */
3189             fix3             = _mm_macc_pd(dx32,fscal,fix3);
3190             fiy3             = _mm_macc_pd(dy32,fscal,fiy3);
3191             fiz3             = _mm_macc_pd(dz32,fscal,fiz3);
3192
3193             fjx2             = _mm_macc_pd(dx32,fscal,fjx2);
3194             fjy2             = _mm_macc_pd(dy32,fscal,fjy2);
3195             fjz2             = _mm_macc_pd(dz32,fscal,fjz2);
3196
3197             }
3198
3199             /**************************
3200              * CALCULATE INTERACTIONS *
3201              **************************/
3202
3203             if (gmx_mm_any_lt(rsq33,rcutoff2))
3204             {
3205
3206             r33              = _mm_mul_pd(rsq33,rinv33);
3207
3208             /* EWALD ELECTROSTATICS */
3209
3210             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
3211             ewrt             = _mm_mul_pd(r33,ewtabscale);
3212             ewitab           = _mm_cvttpd_epi32(ewrt);
3213 #ifdef __XOP__
3214             eweps            = _mm_frcz_pd(ewrt);
3215 #else
3216             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
3217 #endif
3218             twoeweps         = _mm_add_pd(eweps,eweps);
3219             ewitab           = _mm_slli_epi32(ewitab,2);
3220             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
3221             ewtabD           = _mm_setzero_pd();
3222             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
3223             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
3224             ewtabFn          = _mm_setzero_pd();
3225             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
3226             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
3227             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
3228             velec            = _mm_mul_pd(qq33,_mm_sub_pd(rinv33,velec));
3229             felec            = _mm_mul_pd(_mm_mul_pd(qq33,rinv33),_mm_sub_pd(rinvsq33,felec));
3230
3231             d                = _mm_sub_pd(r33,rswitch);
3232             d                = _mm_max_pd(d,_mm_setzero_pd());
3233             d2               = _mm_mul_pd(d,d);
3234             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_macc_pd(d,_mm_macc_pd(d,swV5,swV4),swV3))));
3235
3236             dsw              = _mm_mul_pd(d2,_mm_macc_pd(d,_mm_macc_pd(d,swF4,swF3),swF2));
3237
3238             /* Evaluate switch function */
3239             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
3240             felec            = _mm_msub_pd( felec,sw , _mm_mul_pd(rinv33,_mm_mul_pd(velec,dsw)) );
3241             cutoff_mask      = _mm_cmplt_pd(rsq33,rcutoff2);
3242
3243             fscal            = felec;
3244
3245             fscal            = _mm_and_pd(fscal,cutoff_mask);
3246
3247             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
3248
3249             /* Update vectorial force */
3250             fix3             = _mm_macc_pd(dx33,fscal,fix3);
3251             fiy3             = _mm_macc_pd(dy33,fscal,fiy3);
3252             fiz3             = _mm_macc_pd(dz33,fscal,fiz3);
3253
3254             fjx3             = _mm_macc_pd(dx33,fscal,fjx3);
3255             fjy3             = _mm_macc_pd(dy33,fscal,fjy3);
3256             fjz3             = _mm_macc_pd(dz33,fscal,fjz3);
3257
3258             }
3259
3260             gmx_mm_decrement_4rvec_1ptr_swizzle_pd(f+j_coord_offsetA,fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,fjx2,fjy2,fjz2,fjx3,fjy3,fjz3);
3261
3262             /* Inner loop uses 647 flops */
3263         }
3264
3265         /* End of innermost loop */
3266
3267         gmx_mm_update_iforce_4atom_swizzle_pd(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
3268                                               f+i_coord_offset,fshift+i_shift_offset);
3269
3270         /* Increment number of inner iterations */
3271         inneriter                  += j_index_end - j_index_start;
3272
3273         /* Outer loop uses 24 flops */
3274     }
3275
3276     /* Increment number of outer iterations */
3277     outeriter        += nri;
3278
3279     /* Update outer/inner flops */
3280
3281     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4W4_F,outeriter*24 + inneriter*647);
3282 }