src/gromacs/gmxlib/nonbonded/nb_kernel_avx_128_fma_double/nb_kernel_ElecEwSw_VdwLJSw_GeomW3W3_avx_128_fma_double.c

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