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

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*
  36  * Note: this file was generated by the GROMACS sse2_double kernel generator.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <math.h>
  43
  44 #include "../nb_kernel.h"
  45 #include "types/simple.h"
  46 #include "gromacs/math/vec.h"
  47 #include "nrnb.h"
  48
  49 #include "gromacs/simd/math_x86_sse2_double.h"
  50 #include "kernelutil_x86_sse2_double.h"
  51
  52 /*
  53  * Gromacs nonbonded kernel:   nb_kernel_ElecEwSw_VdwLJSw_GeomW3W3_VF_sse2_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_sse2_double
  61                     (t_nblist                    * gmx_restrict       nlist,
  62                      rvec                        * gmx_restrict          xx,
  63                      rvec                        * gmx_restrict          ff,
  64                      t_forcerec                  * gmx_restrict          fr,
  65                      t_mdatoms                   * gmx_restrict     mdatoms,
  66                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
  67                      t_nrnb                      * gmx_restrict        nrnb)
  68 {
  69     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
  70      * just 0 for non-waters.
  71      * Suffixes A,B refer to j loop unrolling done with SSE double precision, e.g. for the two different
  72      * jnr indices corresponding to data put in the four positions in the SIMD register.
  73      */
  74     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
  75     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
  76     int              jnrA,jnrB;
  77     int              j_coord_offsetA,j_coord_offsetB;
  78     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
  79     real             rcutoff_scalar;
  80     real             *shiftvec,*fshift,*x,*f;
  81     __m128d          tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
  82     int              vdwioffset0;
  83     __m128d          ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
  84     int              vdwioffset1;
  85     __m128d          ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
  86     int              vdwioffset2;
  87     __m128d          ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
  88     int              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,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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
 321             ewitab           = _mm_slli_epi32(ewitab,2);
 322             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
 323             ewtabD           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,1) );
 324             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 325             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
 326             ewtabFn          = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,1) +2);
 327             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 328             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
 329             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
 330             velec            = _mm_mul_pd(qq00,_mm_sub_pd(rinv00,velec));
 331             felec            = _mm_mul_pd(_mm_mul_pd(qq00,rinv00),_mm_sub_pd(rinvsq00,felec));
 332
 333             /* LENNARD-JONES DISPERSION/REPULSION */
 334
 335             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
 336             vvdw6            = _mm_mul_pd(c6_00,rinvsix);
 337             vvdw12           = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
 338             vvdw             = _mm_sub_pd( _mm_mul_pd(vvdw12,one_twelfth) , _mm_mul_pd(vvdw6,one_sixth) );
 339             fvdw             = _mm_mul_pd(_mm_sub_pd(vvdw12,vvdw6),rinvsq00);
 340
 341             d                = _mm_sub_pd(r00,rswitch);
 342             d                = _mm_max_pd(d,_mm_setzero_pd());
 343             d2               = _mm_mul_pd(d,d);
 344             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
 345
 346             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
 347
 348             /* Evaluate switch function */
 349             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 350             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv00,_mm_mul_pd(velec,dsw)) );
 351             fvdw             = _mm_sub_pd( _mm_mul_pd(fvdw,sw) , _mm_mul_pd(rinv00,_mm_mul_pd(vvdw,dsw)) );
 352             velec            = _mm_mul_pd(velec,sw);
 353             vvdw             = _mm_mul_pd(vvdw,sw);
 354             cutoff_mask      = _mm_cmplt_pd(rsq00,rcutoff2);
 355
 356             /* Update potential sum for this i atom from the interaction with this j atom. */
 357             velec            = _mm_and_pd(velec,cutoff_mask);
 358             velecsum         = _mm_add_pd(velecsum,velec);
 359             vvdw             = _mm_and_pd(vvdw,cutoff_mask);
 360             vvdwsum          = _mm_add_pd(vvdwsum,vvdw);
 361
 362             fscal            = _mm_add_pd(felec,fvdw);
 363
 364             fscal            = _mm_and_pd(fscal,cutoff_mask);
 365
 366             /* Calculate temporary vectorial force */
 367             tx               = _mm_mul_pd(fscal,dx00);
 368             ty               = _mm_mul_pd(fscal,dy00);
 369             tz               = _mm_mul_pd(fscal,dz00);
 370
 371             /* Update vectorial force */
 372             fix0             = _mm_add_pd(fix0,tx);
 373             fiy0             = _mm_add_pd(fiy0,ty);
 374             fiz0             = _mm_add_pd(fiz0,tz);
 375
 376             fjx0             = _mm_add_pd(fjx0,tx);
 377             fjy0             = _mm_add_pd(fjy0,ty);
 378             fjz0             = _mm_add_pd(fjz0,tz);
 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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
 397             ewitab           = _mm_slli_epi32(ewitab,2);
 398             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
 399             ewtabD           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,1) );
 400             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 401             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
 402             ewtabFn          = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,1) +2);
 403             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 404             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
 405             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
 406             velec            = _mm_mul_pd(qq01,_mm_sub_pd(rinv01,velec));
 407             felec            = _mm_mul_pd(_mm_mul_pd(qq01,rinv01),_mm_sub_pd(rinvsq01,felec));
 408
 409             d                = _mm_sub_pd(r01,rswitch);
 410             d                = _mm_max_pd(d,_mm_setzero_pd());
 411             d2               = _mm_mul_pd(d,d);
 412             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
 413
 414             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
 415
 416             /* Evaluate switch function */
 417             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 418             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv01,_mm_mul_pd(velec,dsw)) );
 419             velec            = _mm_mul_pd(velec,sw);
 420             cutoff_mask      = _mm_cmplt_pd(rsq01,rcutoff2);
 421
 422             /* Update potential sum for this i atom from the interaction with this j atom. */
 423             velec            = _mm_and_pd(velec,cutoff_mask);
 424             velecsum         = _mm_add_pd(velecsum,velec);
 425
 426             fscal            = felec;
 427
 428             fscal            = _mm_and_pd(fscal,cutoff_mask);
 429
 430             /* Calculate temporary vectorial force */
 431             tx               = _mm_mul_pd(fscal,dx01);
 432             ty               = _mm_mul_pd(fscal,dy01);
 433             tz               = _mm_mul_pd(fscal,dz01);
 434
 435             /* Update vectorial force */
 436             fix0             = _mm_add_pd(fix0,tx);
 437             fiy0             = _mm_add_pd(fiy0,ty);
 438             fiz0             = _mm_add_pd(fiz0,tz);
 439
 440             fjx1             = _mm_add_pd(fjx1,tx);
 441             fjy1             = _mm_add_pd(fjy1,ty);
 442             fjz1             = _mm_add_pd(fjz1,tz);
 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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
 461             ewitab           = _mm_slli_epi32(ewitab,2);
 462             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
 463             ewtabD           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,1) );
 464             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 465             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
 466             ewtabFn          = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,1) +2);
 467             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 468             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
 469             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
 470             velec            = _mm_mul_pd(qq02,_mm_sub_pd(rinv02,velec));
 471             felec            = _mm_mul_pd(_mm_mul_pd(qq02,rinv02),_mm_sub_pd(rinvsq02,felec));
 472
 473             d                = _mm_sub_pd(r02,rswitch);
 474             d                = _mm_max_pd(d,_mm_setzero_pd());
 475             d2               = _mm_mul_pd(d,d);
 476             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
 477
 478             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
 479
 480             /* Evaluate switch function */
 481             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 482             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv02,_mm_mul_pd(velec,dsw)) );
 483             velec            = _mm_mul_pd(velec,sw);
 484             cutoff_mask      = _mm_cmplt_pd(rsq02,rcutoff2);
 485
 486             /* Update potential sum for this i atom from the interaction with this j atom. */
 487             velec            = _mm_and_pd(velec,cutoff_mask);
 488             velecsum         = _mm_add_pd(velecsum,velec);
 489
 490             fscal            = felec;
 491
 492             fscal            = _mm_and_pd(fscal,cutoff_mask);
 493
 494             /* Calculate temporary vectorial force */
 495             tx               = _mm_mul_pd(fscal,dx02);
 496             ty               = _mm_mul_pd(fscal,dy02);
 497             tz               = _mm_mul_pd(fscal,dz02);
 498
 499             /* Update vectorial force */
 500             fix0             = _mm_add_pd(fix0,tx);
 501             fiy0             = _mm_add_pd(fiy0,ty);
 502             fiz0             = _mm_add_pd(fiz0,tz);
 503
 504             fjx2             = _mm_add_pd(fjx2,tx);
 505             fjy2             = _mm_add_pd(fjy2,ty);
 506             fjz2             = _mm_add_pd(fjz2,tz);
 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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
 525             ewitab           = _mm_slli_epi32(ewitab,2);
 526             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
 527             ewtabD           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,1) );
 528             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 529             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
 530             ewtabFn          = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,1) +2);
 531             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 532             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
 533             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
 534             velec            = _mm_mul_pd(qq10,_mm_sub_pd(rinv10,velec));
 535             felec            = _mm_mul_pd(_mm_mul_pd(qq10,rinv10),_mm_sub_pd(rinvsq10,felec));
 536
 537             d                = _mm_sub_pd(r10,rswitch);
 538             d                = _mm_max_pd(d,_mm_setzero_pd());
 539             d2               = _mm_mul_pd(d,d);
 540             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
 541
 542             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
 543
 544             /* Evaluate switch function */
 545             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 546             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv10,_mm_mul_pd(velec,dsw)) );
 547             velec            = _mm_mul_pd(velec,sw);
 548             cutoff_mask      = _mm_cmplt_pd(rsq10,rcutoff2);
 549
 550             /* Update potential sum for this i atom from the interaction with this j atom. */
 551             velec            = _mm_and_pd(velec,cutoff_mask);
 552             velecsum         = _mm_add_pd(velecsum,velec);
 553
 554             fscal            = felec;
 555
 556             fscal            = _mm_and_pd(fscal,cutoff_mask);
 557
 558             /* Calculate temporary vectorial force */
 559             tx               = _mm_mul_pd(fscal,dx10);
 560             ty               = _mm_mul_pd(fscal,dy10);
 561             tz               = _mm_mul_pd(fscal,dz10);
 562
 563             /* Update vectorial force */
 564             fix1             = _mm_add_pd(fix1,tx);
 565             fiy1             = _mm_add_pd(fiy1,ty);
 566             fiz1             = _mm_add_pd(fiz1,tz);
 567
 568             fjx0             = _mm_add_pd(fjx0,tx);
 569             fjy0             = _mm_add_pd(fjy0,ty);
 570             fjz0             = _mm_add_pd(fjz0,tz);
 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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
 589             ewitab           = _mm_slli_epi32(ewitab,2);
 590             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
 591             ewtabD           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,1) );
 592             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 593             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
 594             ewtabFn          = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,1) +2);
 595             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 596             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
 597             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
 598             velec            = _mm_mul_pd(qq11,_mm_sub_pd(rinv11,velec));
 599             felec            = _mm_mul_pd(_mm_mul_pd(qq11,rinv11),_mm_sub_pd(rinvsq11,felec));
 600
 601             d                = _mm_sub_pd(r11,rswitch);
 602             d                = _mm_max_pd(d,_mm_setzero_pd());
 603             d2               = _mm_mul_pd(d,d);
 604             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
 605
 606             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
 607
 608             /* Evaluate switch function */
 609             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 610             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv11,_mm_mul_pd(velec,dsw)) );
 611             velec            = _mm_mul_pd(velec,sw);
 612             cutoff_mask      = _mm_cmplt_pd(rsq11,rcutoff2);
 613
 614             /* Update potential sum for this i atom from the interaction with this j atom. */
 615             velec            = _mm_and_pd(velec,cutoff_mask);
 616             velecsum         = _mm_add_pd(velecsum,velec);
 617
 618             fscal            = felec;
 619
 620             fscal            = _mm_and_pd(fscal,cutoff_mask);
 621
 622             /* Calculate temporary vectorial force */
 623             tx               = _mm_mul_pd(fscal,dx11);
 624             ty               = _mm_mul_pd(fscal,dy11);
 625             tz               = _mm_mul_pd(fscal,dz11);
 626
 627             /* Update vectorial force */
 628             fix1             = _mm_add_pd(fix1,tx);
 629             fiy1             = _mm_add_pd(fiy1,ty);
 630             fiz1             = _mm_add_pd(fiz1,tz);
 631
 632             fjx1             = _mm_add_pd(fjx1,tx);
 633             fjy1             = _mm_add_pd(fjy1,ty);
 634             fjz1             = _mm_add_pd(fjz1,tz);
 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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
 653             ewitab           = _mm_slli_epi32(ewitab,2);
 654             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
 655             ewtabD           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,1) );
 656             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 657             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
 658             ewtabFn          = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,1) +2);
 659             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 660             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
 661             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
 662             velec            = _mm_mul_pd(qq12,_mm_sub_pd(rinv12,velec));
 663             felec            = _mm_mul_pd(_mm_mul_pd(qq12,rinv12),_mm_sub_pd(rinvsq12,felec));
 664
 665             d                = _mm_sub_pd(r12,rswitch);
 666             d                = _mm_max_pd(d,_mm_setzero_pd());
 667             d2               = _mm_mul_pd(d,d);
 668             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
 669
 670             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
 671
 672             /* Evaluate switch function */
 673             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 674             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv12,_mm_mul_pd(velec,dsw)) );
 675             velec            = _mm_mul_pd(velec,sw);
 676             cutoff_mask      = _mm_cmplt_pd(rsq12,rcutoff2);
 677
 678             /* Update potential sum for this i atom from the interaction with this j atom. */
 679             velec            = _mm_and_pd(velec,cutoff_mask);
 680             velecsum         = _mm_add_pd(velecsum,velec);
 681
 682             fscal            = felec;
 683
 684             fscal            = _mm_and_pd(fscal,cutoff_mask);
 685
 686             /* Calculate temporary vectorial force */
 687             tx               = _mm_mul_pd(fscal,dx12);
 688             ty               = _mm_mul_pd(fscal,dy12);
 689             tz               = _mm_mul_pd(fscal,dz12);
 690
 691             /* Update vectorial force */
 692             fix1             = _mm_add_pd(fix1,tx);
 693             fiy1             = _mm_add_pd(fiy1,ty);
 694             fiz1             = _mm_add_pd(fiz1,tz);
 695
 696             fjx2             = _mm_add_pd(fjx2,tx);
 697             fjy2             = _mm_add_pd(fjy2,ty);
 698             fjz2             = _mm_add_pd(fjz2,tz);
 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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
 717             ewitab           = _mm_slli_epi32(ewitab,2);
 718             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
 719             ewtabD           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,1) );
 720             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 721             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
 722             ewtabFn          = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,1) +2);
 723             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 724             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
 725             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
 726             velec            = _mm_mul_pd(qq20,_mm_sub_pd(rinv20,velec));
 727             felec            = _mm_mul_pd(_mm_mul_pd(qq20,rinv20),_mm_sub_pd(rinvsq20,felec));
 728
 729             d                = _mm_sub_pd(r20,rswitch);
 730             d                = _mm_max_pd(d,_mm_setzero_pd());
 731             d2               = _mm_mul_pd(d,d);
 732             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
 733
 734             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
 735
 736             /* Evaluate switch function */
 737             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 738             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv20,_mm_mul_pd(velec,dsw)) );
 739             velec            = _mm_mul_pd(velec,sw);
 740             cutoff_mask      = _mm_cmplt_pd(rsq20,rcutoff2);
 741
 742             /* Update potential sum for this i atom from the interaction with this j atom. */
 743             velec            = _mm_and_pd(velec,cutoff_mask);
 744             velecsum         = _mm_add_pd(velecsum,velec);
 745
 746             fscal            = felec;
 747
 748             fscal            = _mm_and_pd(fscal,cutoff_mask);
 749
 750             /* Calculate temporary vectorial force */
 751             tx               = _mm_mul_pd(fscal,dx20);
 752             ty               = _mm_mul_pd(fscal,dy20);
 753             tz               = _mm_mul_pd(fscal,dz20);
 754
 755             /* Update vectorial force */
 756             fix2             = _mm_add_pd(fix2,tx);
 757             fiy2             = _mm_add_pd(fiy2,ty);
 758             fiz2             = _mm_add_pd(fiz2,tz);
 759
 760             fjx0             = _mm_add_pd(fjx0,tx);
 761             fjy0             = _mm_add_pd(fjy0,ty);
 762             fjz0             = _mm_add_pd(fjz0,tz);
 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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
 781             ewitab           = _mm_slli_epi32(ewitab,2);
 782             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
 783             ewtabD           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,1) );
 784             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 785             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
 786             ewtabFn          = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,1) +2);
 787             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 788             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
 789             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
 790             velec            = _mm_mul_pd(qq21,_mm_sub_pd(rinv21,velec));
 791             felec            = _mm_mul_pd(_mm_mul_pd(qq21,rinv21),_mm_sub_pd(rinvsq21,felec));
 792
 793             d                = _mm_sub_pd(r21,rswitch);
 794             d                = _mm_max_pd(d,_mm_setzero_pd());
 795             d2               = _mm_mul_pd(d,d);
 796             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
 797
 798             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
 799
 800             /* Evaluate switch function */
 801             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 802             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv21,_mm_mul_pd(velec,dsw)) );
 803             velec            = _mm_mul_pd(velec,sw);
 804             cutoff_mask      = _mm_cmplt_pd(rsq21,rcutoff2);
 805
 806             /* Update potential sum for this i atom from the interaction with this j atom. */
 807             velec            = _mm_and_pd(velec,cutoff_mask);
 808             velecsum         = _mm_add_pd(velecsum,velec);
 809
 810             fscal            = felec;
 811
 812             fscal            = _mm_and_pd(fscal,cutoff_mask);
 813
 814             /* Calculate temporary vectorial force */
 815             tx               = _mm_mul_pd(fscal,dx21);
 816             ty               = _mm_mul_pd(fscal,dy21);
 817             tz               = _mm_mul_pd(fscal,dz21);
 818
 819             /* Update vectorial force */
 820             fix2             = _mm_add_pd(fix2,tx);
 821             fiy2             = _mm_add_pd(fiy2,ty);
 822             fiz2             = _mm_add_pd(fiz2,tz);
 823
 824             fjx1             = _mm_add_pd(fjx1,tx);
 825             fjy1             = _mm_add_pd(fjy1,ty);
 826             fjz1             = _mm_add_pd(fjz1,tz);
 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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
 845             ewitab           = _mm_slli_epi32(ewitab,2);
 846             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
 847             ewtabD           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,1) );
 848             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 849             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
 850             ewtabFn          = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,1) +2);
 851             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 852             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
 853             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
 854             velec            = _mm_mul_pd(qq22,_mm_sub_pd(rinv22,velec));
 855             felec            = _mm_mul_pd(_mm_mul_pd(qq22,rinv22),_mm_sub_pd(rinvsq22,felec));
 856
 857             d                = _mm_sub_pd(r22,rswitch);
 858             d                = _mm_max_pd(d,_mm_setzero_pd());
 859             d2               = _mm_mul_pd(d,d);
 860             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
 861
 862             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
 863
 864             /* Evaluate switch function */
 865             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 866             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv22,_mm_mul_pd(velec,dsw)) );
 867             velec            = _mm_mul_pd(velec,sw);
 868             cutoff_mask      = _mm_cmplt_pd(rsq22,rcutoff2);
 869
 870             /* Update potential sum for this i atom from the interaction with this j atom. */
 871             velec            = _mm_and_pd(velec,cutoff_mask);
 872             velecsum         = _mm_add_pd(velecsum,velec);
 873
 874             fscal            = felec;
 875
 876             fscal            = _mm_and_pd(fscal,cutoff_mask);
 877
 878             /* Calculate temporary vectorial force */
 879             tx               = _mm_mul_pd(fscal,dx22);
 880             ty               = _mm_mul_pd(fscal,dy22);
 881             tz               = _mm_mul_pd(fscal,dz22);
 882
 883             /* Update vectorial force */
 884             fix2             = _mm_add_pd(fix2,tx);
 885             fiy2             = _mm_add_pd(fiy2,ty);
 886             fiz2             = _mm_add_pd(fiz2,tz);
 887
 888             fjx2             = _mm_add_pd(fjx2,tx);
 889             fjy2             = _mm_add_pd(fjy2,ty);
 890             fjz2             = _mm_add_pd(fjz2,tz);
 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 603 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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
 994             ewitab           = _mm_slli_epi32(ewitab,2);
 995             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
 996             ewtabD           = _mm_setzero_pd();
 997             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 998             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
 999             ewtabFn          = _mm_setzero_pd();
1000             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1001             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
1002             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
1003             velec            = _mm_mul_pd(qq00,_mm_sub_pd(rinv00,velec));
1004             felec            = _mm_mul_pd(_mm_mul_pd(qq00,rinv00),_mm_sub_pd(rinvsq00,felec));
1005
1006             /* LENNARD-JONES DISPERSION/REPULSION */
1007
1008             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
1009             vvdw6            = _mm_mul_pd(c6_00,rinvsix);
1010             vvdw12           = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
1011             vvdw             = _mm_sub_pd( _mm_mul_pd(vvdw12,one_twelfth) , _mm_mul_pd(vvdw6,one_sixth) );
1012             fvdw             = _mm_mul_pd(_mm_sub_pd(vvdw12,vvdw6),rinvsq00);
1013
1014             d                = _mm_sub_pd(r00,rswitch);
1015             d                = _mm_max_pd(d,_mm_setzero_pd());
1016             d2               = _mm_mul_pd(d,d);
1017             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
1018
1019             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
1020
1021             /* Evaluate switch function */
1022             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1023             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv00,_mm_mul_pd(velec,dsw)) );
1024             fvdw             = _mm_sub_pd( _mm_mul_pd(fvdw,sw) , _mm_mul_pd(rinv00,_mm_mul_pd(vvdw,dsw)) );
1025             velec            = _mm_mul_pd(velec,sw);
1026             vvdw             = _mm_mul_pd(vvdw,sw);
1027             cutoff_mask      = _mm_cmplt_pd(rsq00,rcutoff2);
1028
1029             /* Update potential sum for this i atom from the interaction with this j atom. */
1030             velec            = _mm_and_pd(velec,cutoff_mask);
1031             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
1032             velecsum         = _mm_add_pd(velecsum,velec);
1033             vvdw             = _mm_and_pd(vvdw,cutoff_mask);
1034             vvdw             = _mm_unpacklo_pd(vvdw,_mm_setzero_pd());
1035             vvdwsum          = _mm_add_pd(vvdwsum,vvdw);
1036
1037             fscal            = _mm_add_pd(felec,fvdw);
1038
1039             fscal            = _mm_and_pd(fscal,cutoff_mask);
1040
1041             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1042
1043             /* Calculate temporary vectorial force */
1044             tx               = _mm_mul_pd(fscal,dx00);
1045             ty               = _mm_mul_pd(fscal,dy00);
1046             tz               = _mm_mul_pd(fscal,dz00);
1047
1048             /* Update vectorial force */
1049             fix0             = _mm_add_pd(fix0,tx);
1050             fiy0             = _mm_add_pd(fiy0,ty);
1051             fiz0             = _mm_add_pd(fiz0,tz);
1052
1053             fjx0             = _mm_add_pd(fjx0,tx);
1054             fjy0             = _mm_add_pd(fjy0,ty);
1055             fjz0             = _mm_add_pd(fjz0,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
1074             ewitab           = _mm_slli_epi32(ewitab,2);
1075             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
1076             ewtabD           = _mm_setzero_pd();
1077             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1078             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
1079             ewtabFn          = _mm_setzero_pd();
1080             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1081             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
1082             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
1083             velec            = _mm_mul_pd(qq01,_mm_sub_pd(rinv01,velec));
1084             felec            = _mm_mul_pd(_mm_mul_pd(qq01,rinv01),_mm_sub_pd(rinvsq01,felec));
1085
1086             d                = _mm_sub_pd(r01,rswitch);
1087             d                = _mm_max_pd(d,_mm_setzero_pd());
1088             d2               = _mm_mul_pd(d,d);
1089             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
1090
1091             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
1092
1093             /* Evaluate switch function */
1094             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1095             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv01,_mm_mul_pd(velec,dsw)) );
1096             velec            = _mm_mul_pd(velec,sw);
1097             cutoff_mask      = _mm_cmplt_pd(rsq01,rcutoff2);
1098
1099             /* Update potential sum for this i atom from the interaction with this j atom. */
1100             velec            = _mm_and_pd(velec,cutoff_mask);
1101             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
1102             velecsum         = _mm_add_pd(velecsum,velec);
1103
1104             fscal            = felec;
1105
1106             fscal            = _mm_and_pd(fscal,cutoff_mask);
1107
1108             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1109
1110             /* Calculate temporary vectorial force */
1111             tx               = _mm_mul_pd(fscal,dx01);
1112             ty               = _mm_mul_pd(fscal,dy01);
1113             tz               = _mm_mul_pd(fscal,dz01);
1114
1115             /* Update vectorial force */
1116             fix0             = _mm_add_pd(fix0,tx);
1117             fiy0             = _mm_add_pd(fiy0,ty);
1118             fiz0             = _mm_add_pd(fiz0,tz);
1119
1120             fjx1             = _mm_add_pd(fjx1,tx);
1121             fjy1             = _mm_add_pd(fjy1,ty);
1122             fjz1             = _mm_add_pd(fjz1,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
1141             ewitab           = _mm_slli_epi32(ewitab,2);
1142             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
1143             ewtabD           = _mm_setzero_pd();
1144             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1145             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
1146             ewtabFn          = _mm_setzero_pd();
1147             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1148             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
1149             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
1150             velec            = _mm_mul_pd(qq02,_mm_sub_pd(rinv02,velec));
1151             felec            = _mm_mul_pd(_mm_mul_pd(qq02,rinv02),_mm_sub_pd(rinvsq02,felec));
1152
1153             d                = _mm_sub_pd(r02,rswitch);
1154             d                = _mm_max_pd(d,_mm_setzero_pd());
1155             d2               = _mm_mul_pd(d,d);
1156             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
1157
1158             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
1159
1160             /* Evaluate switch function */
1161             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1162             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv02,_mm_mul_pd(velec,dsw)) );
1163             velec            = _mm_mul_pd(velec,sw);
1164             cutoff_mask      = _mm_cmplt_pd(rsq02,rcutoff2);
1165
1166             /* Update potential sum for this i atom from the interaction with this j atom. */
1167             velec            = _mm_and_pd(velec,cutoff_mask);
1168             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
1169             velecsum         = _mm_add_pd(velecsum,velec);
1170
1171             fscal            = felec;
1172
1173             fscal            = _mm_and_pd(fscal,cutoff_mask);
1174
1175             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1176
1177             /* Calculate temporary vectorial force */
1178             tx               = _mm_mul_pd(fscal,dx02);
1179             ty               = _mm_mul_pd(fscal,dy02);
1180             tz               = _mm_mul_pd(fscal,dz02);
1181
1182             /* Update vectorial force */
1183             fix0             = _mm_add_pd(fix0,tx);
1184             fiy0             = _mm_add_pd(fiy0,ty);
1185             fiz0             = _mm_add_pd(fiz0,tz);
1186
1187             fjx2             = _mm_add_pd(fjx2,tx);
1188             fjy2             = _mm_add_pd(fjy2,ty);
1189             fjz2             = _mm_add_pd(fjz2,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
1208             ewitab           = _mm_slli_epi32(ewitab,2);
1209             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
1210             ewtabD           = _mm_setzero_pd();
1211             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1212             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
1213             ewtabFn          = _mm_setzero_pd();
1214             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1215             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
1216             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
1217             velec            = _mm_mul_pd(qq10,_mm_sub_pd(rinv10,velec));
1218             felec            = _mm_mul_pd(_mm_mul_pd(qq10,rinv10),_mm_sub_pd(rinvsq10,felec));
1219
1220             d                = _mm_sub_pd(r10,rswitch);
1221             d                = _mm_max_pd(d,_mm_setzero_pd());
1222             d2               = _mm_mul_pd(d,d);
1223             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
1224
1225             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
1226
1227             /* Evaluate switch function */
1228             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1229             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv10,_mm_mul_pd(velec,dsw)) );
1230             velec            = _mm_mul_pd(velec,sw);
1231             cutoff_mask      = _mm_cmplt_pd(rsq10,rcutoff2);
1232
1233             /* Update potential sum for this i atom from the interaction with this j atom. */
1234             velec            = _mm_and_pd(velec,cutoff_mask);
1235             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
1236             velecsum         = _mm_add_pd(velecsum,velec);
1237
1238             fscal            = felec;
1239
1240             fscal            = _mm_and_pd(fscal,cutoff_mask);
1241
1242             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1243
1244             /* Calculate temporary vectorial force */
1245             tx               = _mm_mul_pd(fscal,dx10);
1246             ty               = _mm_mul_pd(fscal,dy10);
1247             tz               = _mm_mul_pd(fscal,dz10);
1248
1249             /* Update vectorial force */
1250             fix1             = _mm_add_pd(fix1,tx);
1251             fiy1             = _mm_add_pd(fiy1,ty);
1252             fiz1             = _mm_add_pd(fiz1,tz);
1253
1254             fjx0             = _mm_add_pd(fjx0,tx);
1255             fjy0             = _mm_add_pd(fjy0,ty);
1256             fjz0             = _mm_add_pd(fjz0,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
1275             ewitab           = _mm_slli_epi32(ewitab,2);
1276             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
1277             ewtabD           = _mm_setzero_pd();
1278             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1279             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
1280             ewtabFn          = _mm_setzero_pd();
1281             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1282             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
1283             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
1284             velec            = _mm_mul_pd(qq11,_mm_sub_pd(rinv11,velec));
1285             felec            = _mm_mul_pd(_mm_mul_pd(qq11,rinv11),_mm_sub_pd(rinvsq11,felec));
1286
1287             d                = _mm_sub_pd(r11,rswitch);
1288             d                = _mm_max_pd(d,_mm_setzero_pd());
1289             d2               = _mm_mul_pd(d,d);
1290             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
1291
1292             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
1293
1294             /* Evaluate switch function */
1295             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1296             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv11,_mm_mul_pd(velec,dsw)) );
1297             velec            = _mm_mul_pd(velec,sw);
1298             cutoff_mask      = _mm_cmplt_pd(rsq11,rcutoff2);
1299
1300             /* Update potential sum for this i atom from the interaction with this j atom. */
1301             velec            = _mm_and_pd(velec,cutoff_mask);
1302             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
1303             velecsum         = _mm_add_pd(velecsum,velec);
1304
1305             fscal            = felec;
1306
1307             fscal            = _mm_and_pd(fscal,cutoff_mask);
1308
1309             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1310
1311             /* Calculate temporary vectorial force */
1312             tx               = _mm_mul_pd(fscal,dx11);
1313             ty               = _mm_mul_pd(fscal,dy11);
1314             tz               = _mm_mul_pd(fscal,dz11);
1315
1316             /* Update vectorial force */
1317             fix1             = _mm_add_pd(fix1,tx);
1318             fiy1             = _mm_add_pd(fiy1,ty);
1319             fiz1             = _mm_add_pd(fiz1,tz);
1320
1321             fjx1             = _mm_add_pd(fjx1,tx);
1322             fjy1             = _mm_add_pd(fjy1,ty);
1323             fjz1             = _mm_add_pd(fjz1,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
1342             ewitab           = _mm_slli_epi32(ewitab,2);
1343             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
1344             ewtabD           = _mm_setzero_pd();
1345             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1346             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
1347             ewtabFn          = _mm_setzero_pd();
1348             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1349             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
1350             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
1351             velec            = _mm_mul_pd(qq12,_mm_sub_pd(rinv12,velec));
1352             felec            = _mm_mul_pd(_mm_mul_pd(qq12,rinv12),_mm_sub_pd(rinvsq12,felec));
1353
1354             d                = _mm_sub_pd(r12,rswitch);
1355             d                = _mm_max_pd(d,_mm_setzero_pd());
1356             d2               = _mm_mul_pd(d,d);
1357             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
1358
1359             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
1360
1361             /* Evaluate switch function */
1362             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1363             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv12,_mm_mul_pd(velec,dsw)) );
1364             velec            = _mm_mul_pd(velec,sw);
1365             cutoff_mask      = _mm_cmplt_pd(rsq12,rcutoff2);
1366
1367             /* Update potential sum for this i atom from the interaction with this j atom. */
1368             velec            = _mm_and_pd(velec,cutoff_mask);
1369             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
1370             velecsum         = _mm_add_pd(velecsum,velec);
1371
1372             fscal            = felec;
1373
1374             fscal            = _mm_and_pd(fscal,cutoff_mask);
1375
1376             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1377
1378             /* Calculate temporary vectorial force */
1379             tx               = _mm_mul_pd(fscal,dx12);
1380             ty               = _mm_mul_pd(fscal,dy12);
1381             tz               = _mm_mul_pd(fscal,dz12);
1382
1383             /* Update vectorial force */
1384             fix1             = _mm_add_pd(fix1,tx);
1385             fiy1             = _mm_add_pd(fiy1,ty);
1386             fiz1             = _mm_add_pd(fiz1,tz);
1387
1388             fjx2             = _mm_add_pd(fjx2,tx);
1389             fjy2             = _mm_add_pd(fjy2,ty);
1390             fjz2             = _mm_add_pd(fjz2,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
1409             ewitab           = _mm_slli_epi32(ewitab,2);
1410             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
1411             ewtabD           = _mm_setzero_pd();
1412             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1413             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
1414             ewtabFn          = _mm_setzero_pd();
1415             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1416             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
1417             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
1418             velec            = _mm_mul_pd(qq20,_mm_sub_pd(rinv20,velec));
1419             felec            = _mm_mul_pd(_mm_mul_pd(qq20,rinv20),_mm_sub_pd(rinvsq20,felec));
1420
1421             d                = _mm_sub_pd(r20,rswitch);
1422             d                = _mm_max_pd(d,_mm_setzero_pd());
1423             d2               = _mm_mul_pd(d,d);
1424             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
1425
1426             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
1427
1428             /* Evaluate switch function */
1429             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1430             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv20,_mm_mul_pd(velec,dsw)) );
1431             velec            = _mm_mul_pd(velec,sw);
1432             cutoff_mask      = _mm_cmplt_pd(rsq20,rcutoff2);
1433
1434             /* Update potential sum for this i atom from the interaction with this j atom. */
1435             velec            = _mm_and_pd(velec,cutoff_mask);
1436             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
1437             velecsum         = _mm_add_pd(velecsum,velec);
1438
1439             fscal            = felec;
1440
1441             fscal            = _mm_and_pd(fscal,cutoff_mask);
1442
1443             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1444
1445             /* Calculate temporary vectorial force */
1446             tx               = _mm_mul_pd(fscal,dx20);
1447             ty               = _mm_mul_pd(fscal,dy20);
1448             tz               = _mm_mul_pd(fscal,dz20);
1449
1450             /* Update vectorial force */
1451             fix2             = _mm_add_pd(fix2,tx);
1452             fiy2             = _mm_add_pd(fiy2,ty);
1453             fiz2             = _mm_add_pd(fiz2,tz);
1454
1455             fjx0             = _mm_add_pd(fjx0,tx);
1456             fjy0             = _mm_add_pd(fjy0,ty);
1457             fjz0             = _mm_add_pd(fjz0,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
1476             ewitab           = _mm_slli_epi32(ewitab,2);
1477             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
1478             ewtabD           = _mm_setzero_pd();
1479             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1480             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
1481             ewtabFn          = _mm_setzero_pd();
1482             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1483             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
1484             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
1485             velec            = _mm_mul_pd(qq21,_mm_sub_pd(rinv21,velec));
1486             felec            = _mm_mul_pd(_mm_mul_pd(qq21,rinv21),_mm_sub_pd(rinvsq21,felec));
1487
1488             d                = _mm_sub_pd(r21,rswitch);
1489             d                = _mm_max_pd(d,_mm_setzero_pd());
1490             d2               = _mm_mul_pd(d,d);
1491             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
1492
1493             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
1494
1495             /* Evaluate switch function */
1496             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1497             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv21,_mm_mul_pd(velec,dsw)) );
1498             velec            = _mm_mul_pd(velec,sw);
1499             cutoff_mask      = _mm_cmplt_pd(rsq21,rcutoff2);
1500
1501             /* Update potential sum for this i atom from the interaction with this j atom. */
1502             velec            = _mm_and_pd(velec,cutoff_mask);
1503             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
1504             velecsum         = _mm_add_pd(velecsum,velec);
1505
1506             fscal            = felec;
1507
1508             fscal            = _mm_and_pd(fscal,cutoff_mask);
1509
1510             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1511
1512             /* Calculate temporary vectorial force */
1513             tx               = _mm_mul_pd(fscal,dx21);
1514             ty               = _mm_mul_pd(fscal,dy21);
1515             tz               = _mm_mul_pd(fscal,dz21);
1516
1517             /* Update vectorial force */
1518             fix2             = _mm_add_pd(fix2,tx);
1519             fiy2             = _mm_add_pd(fiy2,ty);
1520             fiz2             = _mm_add_pd(fiz2,tz);
1521
1522             fjx1             = _mm_add_pd(fjx1,tx);
1523             fjy1             = _mm_add_pd(fjy1,ty);
1524             fjz1             = _mm_add_pd(fjz1,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
1543             ewitab           = _mm_slli_epi32(ewitab,2);
1544             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
1545             ewtabD           = _mm_setzero_pd();
1546             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1547             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
1548             ewtabFn          = _mm_setzero_pd();
1549             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1550             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
1551             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
1552             velec            = _mm_mul_pd(qq22,_mm_sub_pd(rinv22,velec));
1553             felec            = _mm_mul_pd(_mm_mul_pd(qq22,rinv22),_mm_sub_pd(rinvsq22,felec));
1554
1555             d                = _mm_sub_pd(r22,rswitch);
1556             d                = _mm_max_pd(d,_mm_setzero_pd());
1557             d2               = _mm_mul_pd(d,d);
1558             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
1559
1560             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
1561
1562             /* Evaluate switch function */
1563             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1564             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv22,_mm_mul_pd(velec,dsw)) );
1565             velec            = _mm_mul_pd(velec,sw);
1566             cutoff_mask      = _mm_cmplt_pd(rsq22,rcutoff2);
1567
1568             /* Update potential sum for this i atom from the interaction with this j atom. */
1569             velec            = _mm_and_pd(velec,cutoff_mask);
1570             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
1571             velecsum         = _mm_add_pd(velecsum,velec);
1572
1573             fscal            = felec;
1574
1575             fscal            = _mm_and_pd(fscal,cutoff_mask);
1576
1577             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1578
1579             /* Calculate temporary vectorial force */
1580             tx               = _mm_mul_pd(fscal,dx22);
1581             ty               = _mm_mul_pd(fscal,dy22);
1582             tz               = _mm_mul_pd(fscal,dz22);
1583
1584             /* Update vectorial force */
1585             fix2             = _mm_add_pd(fix2,tx);
1586             fiy2             = _mm_add_pd(fiy2,ty);
1587             fiz2             = _mm_add_pd(fiz2,tz);
1588
1589             fjx2             = _mm_add_pd(fjx2,tx);
1590             fjy2             = _mm_add_pd(fjy2,ty);
1591             fjz2             = _mm_add_pd(fjz2,tz);
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 603 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*603);
1622 }
1623 /*
1624  * Gromacs nonbonded kernel:   nb_kernel_ElecEwSw_VdwLJSw_GeomW3W3_F_sse2_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_sse2_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,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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
1888             ewitab           = _mm_slli_epi32(ewitab,2);
1889             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
1890             ewtabD           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,1) );
1891             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1892             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
1893             ewtabFn          = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,1) +2);
1894             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1895             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
1896             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
1897             velec            = _mm_mul_pd(qq00,_mm_sub_pd(rinv00,velec));
1898             felec            = _mm_mul_pd(_mm_mul_pd(qq00,rinv00),_mm_sub_pd(rinvsq00,felec));
1899
1900             /* LENNARD-JONES DISPERSION/REPULSION */
1901
1902             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
1903             vvdw6            = _mm_mul_pd(c6_00,rinvsix);
1904             vvdw12           = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
1905             vvdw             = _mm_sub_pd( _mm_mul_pd(vvdw12,one_twelfth) , _mm_mul_pd(vvdw6,one_sixth) );
1906             fvdw             = _mm_mul_pd(_mm_sub_pd(vvdw12,vvdw6),rinvsq00);
1907
1908             d                = _mm_sub_pd(r00,rswitch);
1909             d                = _mm_max_pd(d,_mm_setzero_pd());
1910             d2               = _mm_mul_pd(d,d);
1911             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
1912
1913             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
1914
1915             /* Evaluate switch function */
1916             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1917             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv00,_mm_mul_pd(velec,dsw)) );
1918             fvdw             = _mm_sub_pd( _mm_mul_pd(fvdw,sw) , _mm_mul_pd(rinv00,_mm_mul_pd(vvdw,dsw)) );
1919             cutoff_mask      = _mm_cmplt_pd(rsq00,rcutoff2);
1920
1921             fscal            = _mm_add_pd(felec,fvdw);
1922
1923             fscal            = _mm_and_pd(fscal,cutoff_mask);
1924
1925             /* Calculate temporary vectorial force */
1926             tx               = _mm_mul_pd(fscal,dx00);
1927             ty               = _mm_mul_pd(fscal,dy00);
1928             tz               = _mm_mul_pd(fscal,dz00);
1929
1930             /* Update vectorial force */
1931             fix0             = _mm_add_pd(fix0,tx);
1932             fiy0             = _mm_add_pd(fiy0,ty);
1933             fiz0             = _mm_add_pd(fiz0,tz);
1934
1935             fjx0             = _mm_add_pd(fjx0,tx);
1936             fjy0             = _mm_add_pd(fjy0,ty);
1937             fjz0             = _mm_add_pd(fjz0,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
1956             ewitab           = _mm_slli_epi32(ewitab,2);
1957             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
1958             ewtabD           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,1) );
1959             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
1960             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
1961             ewtabFn          = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,1) +2);
1962             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
1963             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
1964             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
1965             velec            = _mm_mul_pd(qq01,_mm_sub_pd(rinv01,velec));
1966             felec            = _mm_mul_pd(_mm_mul_pd(qq01,rinv01),_mm_sub_pd(rinvsq01,felec));
1967
1968             d                = _mm_sub_pd(r01,rswitch);
1969             d                = _mm_max_pd(d,_mm_setzero_pd());
1970             d2               = _mm_mul_pd(d,d);
1971             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
1972
1973             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
1974
1975             /* Evaluate switch function */
1976             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1977             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv01,_mm_mul_pd(velec,dsw)) );
1978             cutoff_mask      = _mm_cmplt_pd(rsq01,rcutoff2);
1979
1980             fscal            = felec;
1981
1982             fscal            = _mm_and_pd(fscal,cutoff_mask);
1983
1984             /* Calculate temporary vectorial force */
1985             tx               = _mm_mul_pd(fscal,dx01);
1986             ty               = _mm_mul_pd(fscal,dy01);
1987             tz               = _mm_mul_pd(fscal,dz01);
1988
1989             /* Update vectorial force */
1990             fix0             = _mm_add_pd(fix0,tx);
1991             fiy0             = _mm_add_pd(fiy0,ty);
1992             fiz0             = _mm_add_pd(fiz0,tz);
1993
1994             fjx1             = _mm_add_pd(fjx1,tx);
1995             fjy1             = _mm_add_pd(fjy1,ty);
1996             fjz1             = _mm_add_pd(fjz1,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
2015             ewitab           = _mm_slli_epi32(ewitab,2);
2016             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2017             ewtabD           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2018             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2019             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
2020             ewtabFn          = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,1) +2);
2021             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2022             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
2023             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
2024             velec            = _mm_mul_pd(qq02,_mm_sub_pd(rinv02,velec));
2025             felec            = _mm_mul_pd(_mm_mul_pd(qq02,rinv02),_mm_sub_pd(rinvsq02,felec));
2026
2027             d                = _mm_sub_pd(r02,rswitch);
2028             d                = _mm_max_pd(d,_mm_setzero_pd());
2029             d2               = _mm_mul_pd(d,d);
2030             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
2031
2032             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
2033
2034             /* Evaluate switch function */
2035             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2036             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv02,_mm_mul_pd(velec,dsw)) );
2037             cutoff_mask      = _mm_cmplt_pd(rsq02,rcutoff2);
2038
2039             fscal            = felec;
2040
2041             fscal            = _mm_and_pd(fscal,cutoff_mask);
2042
2043             /* Calculate temporary vectorial force */
2044             tx               = _mm_mul_pd(fscal,dx02);
2045             ty               = _mm_mul_pd(fscal,dy02);
2046             tz               = _mm_mul_pd(fscal,dz02);
2047
2048             /* Update vectorial force */
2049             fix0             = _mm_add_pd(fix0,tx);
2050             fiy0             = _mm_add_pd(fiy0,ty);
2051             fiz0             = _mm_add_pd(fiz0,tz);
2052
2053             fjx2             = _mm_add_pd(fjx2,tx);
2054             fjy2             = _mm_add_pd(fjy2,ty);
2055             fjz2             = _mm_add_pd(fjz2,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
2074             ewitab           = _mm_slli_epi32(ewitab,2);
2075             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2076             ewtabD           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2077             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2078             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
2079             ewtabFn          = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,1) +2);
2080             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2081             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
2082             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
2083             velec            = _mm_mul_pd(qq10,_mm_sub_pd(rinv10,velec));
2084             felec            = _mm_mul_pd(_mm_mul_pd(qq10,rinv10),_mm_sub_pd(rinvsq10,felec));
2085
2086             d                = _mm_sub_pd(r10,rswitch);
2087             d                = _mm_max_pd(d,_mm_setzero_pd());
2088             d2               = _mm_mul_pd(d,d);
2089             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
2090
2091             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
2092
2093             /* Evaluate switch function */
2094             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2095             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv10,_mm_mul_pd(velec,dsw)) );
2096             cutoff_mask      = _mm_cmplt_pd(rsq10,rcutoff2);
2097
2098             fscal            = felec;
2099
2100             fscal            = _mm_and_pd(fscal,cutoff_mask);
2101
2102             /* Calculate temporary vectorial force */
2103             tx               = _mm_mul_pd(fscal,dx10);
2104             ty               = _mm_mul_pd(fscal,dy10);
2105             tz               = _mm_mul_pd(fscal,dz10);
2106
2107             /* Update vectorial force */
2108             fix1             = _mm_add_pd(fix1,tx);
2109             fiy1             = _mm_add_pd(fiy1,ty);
2110             fiz1             = _mm_add_pd(fiz1,tz);
2111
2112             fjx0             = _mm_add_pd(fjx0,tx);
2113             fjy0             = _mm_add_pd(fjy0,ty);
2114             fjz0             = _mm_add_pd(fjz0,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
2133             ewitab           = _mm_slli_epi32(ewitab,2);
2134             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2135             ewtabD           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2136             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2137             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
2138             ewtabFn          = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,1) +2);
2139             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2140             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
2141             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
2142             velec            = _mm_mul_pd(qq11,_mm_sub_pd(rinv11,velec));
2143             felec            = _mm_mul_pd(_mm_mul_pd(qq11,rinv11),_mm_sub_pd(rinvsq11,felec));
2144
2145             d                = _mm_sub_pd(r11,rswitch);
2146             d                = _mm_max_pd(d,_mm_setzero_pd());
2147             d2               = _mm_mul_pd(d,d);
2148             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
2149
2150             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
2151
2152             /* Evaluate switch function */
2153             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2154             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv11,_mm_mul_pd(velec,dsw)) );
2155             cutoff_mask      = _mm_cmplt_pd(rsq11,rcutoff2);
2156
2157             fscal            = felec;
2158
2159             fscal            = _mm_and_pd(fscal,cutoff_mask);
2160
2161             /* Calculate temporary vectorial force */
2162             tx               = _mm_mul_pd(fscal,dx11);
2163             ty               = _mm_mul_pd(fscal,dy11);
2164             tz               = _mm_mul_pd(fscal,dz11);
2165
2166             /* Update vectorial force */
2167             fix1             = _mm_add_pd(fix1,tx);
2168             fiy1             = _mm_add_pd(fiy1,ty);
2169             fiz1             = _mm_add_pd(fiz1,tz);
2170
2171             fjx1             = _mm_add_pd(fjx1,tx);
2172             fjy1             = _mm_add_pd(fjy1,ty);
2173             fjz1             = _mm_add_pd(fjz1,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
2192             ewitab           = _mm_slli_epi32(ewitab,2);
2193             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2194             ewtabD           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2195             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2196             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
2197             ewtabFn          = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,1) +2);
2198             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2199             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
2200             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
2201             velec            = _mm_mul_pd(qq12,_mm_sub_pd(rinv12,velec));
2202             felec            = _mm_mul_pd(_mm_mul_pd(qq12,rinv12),_mm_sub_pd(rinvsq12,felec));
2203
2204             d                = _mm_sub_pd(r12,rswitch);
2205             d                = _mm_max_pd(d,_mm_setzero_pd());
2206             d2               = _mm_mul_pd(d,d);
2207             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
2208
2209             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
2210
2211             /* Evaluate switch function */
2212             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2213             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv12,_mm_mul_pd(velec,dsw)) );
2214             cutoff_mask      = _mm_cmplt_pd(rsq12,rcutoff2);
2215
2216             fscal            = felec;
2217
2218             fscal            = _mm_and_pd(fscal,cutoff_mask);
2219
2220             /* Calculate temporary vectorial force */
2221             tx               = _mm_mul_pd(fscal,dx12);
2222             ty               = _mm_mul_pd(fscal,dy12);
2223             tz               = _mm_mul_pd(fscal,dz12);
2224
2225             /* Update vectorial force */
2226             fix1             = _mm_add_pd(fix1,tx);
2227             fiy1             = _mm_add_pd(fiy1,ty);
2228             fiz1             = _mm_add_pd(fiz1,tz);
2229
2230             fjx2             = _mm_add_pd(fjx2,tx);
2231             fjy2             = _mm_add_pd(fjy2,ty);
2232             fjz2             = _mm_add_pd(fjz2,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
2251             ewitab           = _mm_slli_epi32(ewitab,2);
2252             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2253             ewtabD           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2254             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2255             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
2256             ewtabFn          = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,1) +2);
2257             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2258             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
2259             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
2260             velec            = _mm_mul_pd(qq20,_mm_sub_pd(rinv20,velec));
2261             felec            = _mm_mul_pd(_mm_mul_pd(qq20,rinv20),_mm_sub_pd(rinvsq20,felec));
2262
2263             d                = _mm_sub_pd(r20,rswitch);
2264             d                = _mm_max_pd(d,_mm_setzero_pd());
2265             d2               = _mm_mul_pd(d,d);
2266             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
2267
2268             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
2269
2270             /* Evaluate switch function */
2271             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2272             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv20,_mm_mul_pd(velec,dsw)) );
2273             cutoff_mask      = _mm_cmplt_pd(rsq20,rcutoff2);
2274
2275             fscal            = felec;
2276
2277             fscal            = _mm_and_pd(fscal,cutoff_mask);
2278
2279             /* Calculate temporary vectorial force */
2280             tx               = _mm_mul_pd(fscal,dx20);
2281             ty               = _mm_mul_pd(fscal,dy20);
2282             tz               = _mm_mul_pd(fscal,dz20);
2283
2284             /* Update vectorial force */
2285             fix2             = _mm_add_pd(fix2,tx);
2286             fiy2             = _mm_add_pd(fiy2,ty);
2287             fiz2             = _mm_add_pd(fiz2,tz);
2288
2289             fjx0             = _mm_add_pd(fjx0,tx);
2290             fjy0             = _mm_add_pd(fjy0,ty);
2291             fjz0             = _mm_add_pd(fjz0,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
2310             ewitab           = _mm_slli_epi32(ewitab,2);
2311             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2312             ewtabD           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2313             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2314             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
2315             ewtabFn          = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,1) +2);
2316             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2317             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
2318             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
2319             velec            = _mm_mul_pd(qq21,_mm_sub_pd(rinv21,velec));
2320             felec            = _mm_mul_pd(_mm_mul_pd(qq21,rinv21),_mm_sub_pd(rinvsq21,felec));
2321
2322             d                = _mm_sub_pd(r21,rswitch);
2323             d                = _mm_max_pd(d,_mm_setzero_pd());
2324             d2               = _mm_mul_pd(d,d);
2325             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
2326
2327             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
2328
2329             /* Evaluate switch function */
2330             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2331             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv21,_mm_mul_pd(velec,dsw)) );
2332             cutoff_mask      = _mm_cmplt_pd(rsq21,rcutoff2);
2333
2334             fscal            = felec;
2335
2336             fscal            = _mm_and_pd(fscal,cutoff_mask);
2337
2338             /* Calculate temporary vectorial force */
2339             tx               = _mm_mul_pd(fscal,dx21);
2340             ty               = _mm_mul_pd(fscal,dy21);
2341             tz               = _mm_mul_pd(fscal,dz21);
2342
2343             /* Update vectorial force */
2344             fix2             = _mm_add_pd(fix2,tx);
2345             fiy2             = _mm_add_pd(fiy2,ty);
2346             fiz2             = _mm_add_pd(fiz2,tz);
2347
2348             fjx1             = _mm_add_pd(fjx1,tx);
2349             fjy1             = _mm_add_pd(fjy1,ty);
2350             fjz1             = _mm_add_pd(fjz1,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
2369             ewitab           = _mm_slli_epi32(ewitab,2);
2370             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2371             ewtabD           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2372             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2373             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
2374             ewtabFn          = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,1) +2);
2375             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2376             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
2377             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
2378             velec            = _mm_mul_pd(qq22,_mm_sub_pd(rinv22,velec));
2379             felec            = _mm_mul_pd(_mm_mul_pd(qq22,rinv22),_mm_sub_pd(rinvsq22,felec));
2380
2381             d                = _mm_sub_pd(r22,rswitch);
2382             d                = _mm_max_pd(d,_mm_setzero_pd());
2383             d2               = _mm_mul_pd(d,d);
2384             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
2385
2386             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
2387
2388             /* Evaluate switch function */
2389             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2390             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv22,_mm_mul_pd(velec,dsw)) );
2391             cutoff_mask      = _mm_cmplt_pd(rsq22,rcutoff2);
2392
2393             fscal            = felec;
2394
2395             fscal            = _mm_and_pd(fscal,cutoff_mask);
2396
2397             /* Calculate temporary vectorial force */
2398             tx               = _mm_mul_pd(fscal,dx22);
2399             ty               = _mm_mul_pd(fscal,dy22);
2400             tz               = _mm_mul_pd(fscal,dz22);
2401
2402             /* Update vectorial force */
2403             fix2             = _mm_add_pd(fix2,tx);
2404             fiy2             = _mm_add_pd(fiy2,ty);
2405             fiz2             = _mm_add_pd(fiz2,tz);
2406
2407             fjx2             = _mm_add_pd(fjx2,tx);
2408             fjy2             = _mm_add_pd(fjy2,ty);
2409             fjz2             = _mm_add_pd(fjz2,tz);
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 573 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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
2513             ewitab           = _mm_slli_epi32(ewitab,2);
2514             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2515             ewtabD           = _mm_setzero_pd();
2516             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2517             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
2518             ewtabFn          = _mm_setzero_pd();
2519             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2520             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
2521             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
2522             velec            = _mm_mul_pd(qq00,_mm_sub_pd(rinv00,velec));
2523             felec            = _mm_mul_pd(_mm_mul_pd(qq00,rinv00),_mm_sub_pd(rinvsq00,felec));
2524
2525             /* LENNARD-JONES DISPERSION/REPULSION */
2526
2527             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
2528             vvdw6            = _mm_mul_pd(c6_00,rinvsix);
2529             vvdw12           = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
2530             vvdw             = _mm_sub_pd( _mm_mul_pd(vvdw12,one_twelfth) , _mm_mul_pd(vvdw6,one_sixth) );
2531             fvdw             = _mm_mul_pd(_mm_sub_pd(vvdw12,vvdw6),rinvsq00);
2532
2533             d                = _mm_sub_pd(r00,rswitch);
2534             d                = _mm_max_pd(d,_mm_setzero_pd());
2535             d2               = _mm_mul_pd(d,d);
2536             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
2537
2538             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
2539
2540             /* Evaluate switch function */
2541             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2542             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv00,_mm_mul_pd(velec,dsw)) );
2543             fvdw             = _mm_sub_pd( _mm_mul_pd(fvdw,sw) , _mm_mul_pd(rinv00,_mm_mul_pd(vvdw,dsw)) );
2544             cutoff_mask      = _mm_cmplt_pd(rsq00,rcutoff2);
2545
2546             fscal            = _mm_add_pd(felec,fvdw);
2547
2548             fscal            = _mm_and_pd(fscal,cutoff_mask);
2549
2550             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
2551
2552             /* Calculate temporary vectorial force */
2553             tx               = _mm_mul_pd(fscal,dx00);
2554             ty               = _mm_mul_pd(fscal,dy00);
2555             tz               = _mm_mul_pd(fscal,dz00);
2556
2557             /* Update vectorial force */
2558             fix0             = _mm_add_pd(fix0,tx);
2559             fiy0             = _mm_add_pd(fiy0,ty);
2560             fiz0             = _mm_add_pd(fiz0,tz);
2561
2562             fjx0             = _mm_add_pd(fjx0,tx);
2563             fjy0             = _mm_add_pd(fjy0,ty);
2564             fjz0             = _mm_add_pd(fjz0,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
2583             ewitab           = _mm_slli_epi32(ewitab,2);
2584             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2585             ewtabD           = _mm_setzero_pd();
2586             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2587             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
2588             ewtabFn          = _mm_setzero_pd();
2589             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2590             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
2591             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
2592             velec            = _mm_mul_pd(qq01,_mm_sub_pd(rinv01,velec));
2593             felec            = _mm_mul_pd(_mm_mul_pd(qq01,rinv01),_mm_sub_pd(rinvsq01,felec));
2594
2595             d                = _mm_sub_pd(r01,rswitch);
2596             d                = _mm_max_pd(d,_mm_setzero_pd());
2597             d2               = _mm_mul_pd(d,d);
2598             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
2599
2600             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
2601
2602             /* Evaluate switch function */
2603             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2604             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv01,_mm_mul_pd(velec,dsw)) );
2605             cutoff_mask      = _mm_cmplt_pd(rsq01,rcutoff2);
2606
2607             fscal            = felec;
2608
2609             fscal            = _mm_and_pd(fscal,cutoff_mask);
2610
2611             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
2612
2613             /* Calculate temporary vectorial force */
2614             tx               = _mm_mul_pd(fscal,dx01);
2615             ty               = _mm_mul_pd(fscal,dy01);
2616             tz               = _mm_mul_pd(fscal,dz01);
2617
2618             /* Update vectorial force */
2619             fix0             = _mm_add_pd(fix0,tx);
2620             fiy0             = _mm_add_pd(fiy0,ty);
2621             fiz0             = _mm_add_pd(fiz0,tz);
2622
2623             fjx1             = _mm_add_pd(fjx1,tx);
2624             fjy1             = _mm_add_pd(fjy1,ty);
2625             fjz1             = _mm_add_pd(fjz1,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
2644             ewitab           = _mm_slli_epi32(ewitab,2);
2645             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2646             ewtabD           = _mm_setzero_pd();
2647             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2648             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
2649             ewtabFn          = _mm_setzero_pd();
2650             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2651             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
2652             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
2653             velec            = _mm_mul_pd(qq02,_mm_sub_pd(rinv02,velec));
2654             felec            = _mm_mul_pd(_mm_mul_pd(qq02,rinv02),_mm_sub_pd(rinvsq02,felec));
2655
2656             d                = _mm_sub_pd(r02,rswitch);
2657             d                = _mm_max_pd(d,_mm_setzero_pd());
2658             d2               = _mm_mul_pd(d,d);
2659             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
2660
2661             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
2662
2663             /* Evaluate switch function */
2664             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2665             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv02,_mm_mul_pd(velec,dsw)) );
2666             cutoff_mask      = _mm_cmplt_pd(rsq02,rcutoff2);
2667
2668             fscal            = felec;
2669
2670             fscal            = _mm_and_pd(fscal,cutoff_mask);
2671
2672             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
2673
2674             /* Calculate temporary vectorial force */
2675             tx               = _mm_mul_pd(fscal,dx02);
2676             ty               = _mm_mul_pd(fscal,dy02);
2677             tz               = _mm_mul_pd(fscal,dz02);
2678
2679             /* Update vectorial force */
2680             fix0             = _mm_add_pd(fix0,tx);
2681             fiy0             = _mm_add_pd(fiy0,ty);
2682             fiz0             = _mm_add_pd(fiz0,tz);
2683
2684             fjx2             = _mm_add_pd(fjx2,tx);
2685             fjy2             = _mm_add_pd(fjy2,ty);
2686             fjz2             = _mm_add_pd(fjz2,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
2705             ewitab           = _mm_slli_epi32(ewitab,2);
2706             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2707             ewtabD           = _mm_setzero_pd();
2708             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2709             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
2710             ewtabFn          = _mm_setzero_pd();
2711             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2712             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
2713             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
2714             velec            = _mm_mul_pd(qq10,_mm_sub_pd(rinv10,velec));
2715             felec            = _mm_mul_pd(_mm_mul_pd(qq10,rinv10),_mm_sub_pd(rinvsq10,felec));
2716
2717             d                = _mm_sub_pd(r10,rswitch);
2718             d                = _mm_max_pd(d,_mm_setzero_pd());
2719             d2               = _mm_mul_pd(d,d);
2720             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
2721
2722             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
2723
2724             /* Evaluate switch function */
2725             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2726             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv10,_mm_mul_pd(velec,dsw)) );
2727             cutoff_mask      = _mm_cmplt_pd(rsq10,rcutoff2);
2728
2729             fscal            = felec;
2730
2731             fscal            = _mm_and_pd(fscal,cutoff_mask);
2732
2733             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
2734
2735             /* Calculate temporary vectorial force */
2736             tx               = _mm_mul_pd(fscal,dx10);
2737             ty               = _mm_mul_pd(fscal,dy10);
2738             tz               = _mm_mul_pd(fscal,dz10);
2739
2740             /* Update vectorial force */
2741             fix1             = _mm_add_pd(fix1,tx);
2742             fiy1             = _mm_add_pd(fiy1,ty);
2743             fiz1             = _mm_add_pd(fiz1,tz);
2744
2745             fjx0             = _mm_add_pd(fjx0,tx);
2746             fjy0             = _mm_add_pd(fjy0,ty);
2747             fjz0             = _mm_add_pd(fjz0,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
2766             ewitab           = _mm_slli_epi32(ewitab,2);
2767             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2768             ewtabD           = _mm_setzero_pd();
2769             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2770             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
2771             ewtabFn          = _mm_setzero_pd();
2772             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2773             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
2774             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
2775             velec            = _mm_mul_pd(qq11,_mm_sub_pd(rinv11,velec));
2776             felec            = _mm_mul_pd(_mm_mul_pd(qq11,rinv11),_mm_sub_pd(rinvsq11,felec));
2777
2778             d                = _mm_sub_pd(r11,rswitch);
2779             d                = _mm_max_pd(d,_mm_setzero_pd());
2780             d2               = _mm_mul_pd(d,d);
2781             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
2782
2783             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
2784
2785             /* Evaluate switch function */
2786             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2787             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv11,_mm_mul_pd(velec,dsw)) );
2788             cutoff_mask      = _mm_cmplt_pd(rsq11,rcutoff2);
2789
2790             fscal            = felec;
2791
2792             fscal            = _mm_and_pd(fscal,cutoff_mask);
2793
2794             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
2795
2796             /* Calculate temporary vectorial force */
2797             tx               = _mm_mul_pd(fscal,dx11);
2798             ty               = _mm_mul_pd(fscal,dy11);
2799             tz               = _mm_mul_pd(fscal,dz11);
2800
2801             /* Update vectorial force */
2802             fix1             = _mm_add_pd(fix1,tx);
2803             fiy1             = _mm_add_pd(fiy1,ty);
2804             fiz1             = _mm_add_pd(fiz1,tz);
2805
2806             fjx1             = _mm_add_pd(fjx1,tx);
2807             fjy1             = _mm_add_pd(fjy1,ty);
2808             fjz1             = _mm_add_pd(fjz1,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
2827             ewitab           = _mm_slli_epi32(ewitab,2);
2828             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2829             ewtabD           = _mm_setzero_pd();
2830             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2831             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
2832             ewtabFn          = _mm_setzero_pd();
2833             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2834             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
2835             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
2836             velec            = _mm_mul_pd(qq12,_mm_sub_pd(rinv12,velec));
2837             felec            = _mm_mul_pd(_mm_mul_pd(qq12,rinv12),_mm_sub_pd(rinvsq12,felec));
2838
2839             d                = _mm_sub_pd(r12,rswitch);
2840             d                = _mm_max_pd(d,_mm_setzero_pd());
2841             d2               = _mm_mul_pd(d,d);
2842             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
2843
2844             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
2845
2846             /* Evaluate switch function */
2847             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2848             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv12,_mm_mul_pd(velec,dsw)) );
2849             cutoff_mask      = _mm_cmplt_pd(rsq12,rcutoff2);
2850
2851             fscal            = felec;
2852
2853             fscal            = _mm_and_pd(fscal,cutoff_mask);
2854
2855             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
2856
2857             /* Calculate temporary vectorial force */
2858             tx               = _mm_mul_pd(fscal,dx12);
2859             ty               = _mm_mul_pd(fscal,dy12);
2860             tz               = _mm_mul_pd(fscal,dz12);
2861
2862             /* Update vectorial force */
2863             fix1             = _mm_add_pd(fix1,tx);
2864             fiy1             = _mm_add_pd(fiy1,ty);
2865             fiz1             = _mm_add_pd(fiz1,tz);
2866
2867             fjx2             = _mm_add_pd(fjx2,tx);
2868             fjy2             = _mm_add_pd(fjy2,ty);
2869             fjz2             = _mm_add_pd(fjz2,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
2888             ewitab           = _mm_slli_epi32(ewitab,2);
2889             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2890             ewtabD           = _mm_setzero_pd();
2891             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2892             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
2893             ewtabFn          = _mm_setzero_pd();
2894             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2895             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
2896             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
2897             velec            = _mm_mul_pd(qq20,_mm_sub_pd(rinv20,velec));
2898             felec            = _mm_mul_pd(_mm_mul_pd(qq20,rinv20),_mm_sub_pd(rinvsq20,felec));
2899
2900             d                = _mm_sub_pd(r20,rswitch);
2901             d                = _mm_max_pd(d,_mm_setzero_pd());
2902             d2               = _mm_mul_pd(d,d);
2903             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
2904
2905             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
2906
2907             /* Evaluate switch function */
2908             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2909             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv20,_mm_mul_pd(velec,dsw)) );
2910             cutoff_mask      = _mm_cmplt_pd(rsq20,rcutoff2);
2911
2912             fscal            = felec;
2913
2914             fscal            = _mm_and_pd(fscal,cutoff_mask);
2915
2916             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
2917
2918             /* Calculate temporary vectorial force */
2919             tx               = _mm_mul_pd(fscal,dx20);
2920             ty               = _mm_mul_pd(fscal,dy20);
2921             tz               = _mm_mul_pd(fscal,dz20);
2922
2923             /* Update vectorial force */
2924             fix2             = _mm_add_pd(fix2,tx);
2925             fiy2             = _mm_add_pd(fiy2,ty);
2926             fiz2             = _mm_add_pd(fiz2,tz);
2927
2928             fjx0             = _mm_add_pd(fjx0,tx);
2929             fjy0             = _mm_add_pd(fjy0,ty);
2930             fjz0             = _mm_add_pd(fjz0,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
2949             ewitab           = _mm_slli_epi32(ewitab,2);
2950             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2951             ewtabD           = _mm_setzero_pd();
2952             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
2953             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
2954             ewtabFn          = _mm_setzero_pd();
2955             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
2956             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
2957             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
2958             velec            = _mm_mul_pd(qq21,_mm_sub_pd(rinv21,velec));
2959             felec            = _mm_mul_pd(_mm_mul_pd(qq21,rinv21),_mm_sub_pd(rinvsq21,felec));
2960
2961             d                = _mm_sub_pd(r21,rswitch);
2962             d                = _mm_max_pd(d,_mm_setzero_pd());
2963             d2               = _mm_mul_pd(d,d);
2964             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
2965
2966             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
2967
2968             /* Evaluate switch function */
2969             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2970             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv21,_mm_mul_pd(velec,dsw)) );
2971             cutoff_mask      = _mm_cmplt_pd(rsq21,rcutoff2);
2972
2973             fscal            = felec;
2974
2975             fscal            = _mm_and_pd(fscal,cutoff_mask);
2976
2977             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
2978
2979             /* Calculate temporary vectorial force */
2980             tx               = _mm_mul_pd(fscal,dx21);
2981             ty               = _mm_mul_pd(fscal,dy21);
2982             tz               = _mm_mul_pd(fscal,dz21);
2983
2984             /* Update vectorial force */
2985             fix2             = _mm_add_pd(fix2,tx);
2986             fiy2             = _mm_add_pd(fiy2,ty);
2987             fiz2             = _mm_add_pd(fiz2,tz);
2988
2989             fjx1             = _mm_add_pd(fjx1,tx);
2990             fjy1             = _mm_add_pd(fjy1,ty);
2991             fjz1             = _mm_add_pd(fjz1,tz);
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             eweps            = _mm_sub_pd(ewrt,_mm_cvtepi32_pd(ewitab));
3010             ewitab           = _mm_slli_epi32(ewitab,2);
3011             ewtabF           = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
3012             ewtabD           = _mm_setzero_pd();
3013             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
3014             ewtabV           = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
3015             ewtabFn          = _mm_setzero_pd();
3016             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
3017             felec            = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
3018             velec            = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
3019             velec            = _mm_mul_pd(qq22,_mm_sub_pd(rinv22,velec));
3020             felec            = _mm_mul_pd(_mm_mul_pd(qq22,rinv22),_mm_sub_pd(rinvsq22,felec));
3021
3022             d                = _mm_sub_pd(r22,rswitch);
3023             d                = _mm_max_pd(d,_mm_setzero_pd());
3024             d2               = _mm_mul_pd(d,d);
3025             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
3026
3027             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
3028
3029             /* Evaluate switch function */
3030             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
3031             felec            = _mm_sub_pd( _mm_mul_pd(felec,sw) , _mm_mul_pd(rinv22,_mm_mul_pd(velec,dsw)) );
3032             cutoff_mask      = _mm_cmplt_pd(rsq22,rcutoff2);
3033
3034             fscal            = felec;
3035
3036             fscal            = _mm_and_pd(fscal,cutoff_mask);
3037
3038             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
3039
3040             /* Calculate temporary vectorial force */
3041             tx               = _mm_mul_pd(fscal,dx22);
3042             ty               = _mm_mul_pd(fscal,dy22);
3043             tz               = _mm_mul_pd(fscal,dz22);
3044
3045             /* Update vectorial force */
3046             fix2             = _mm_add_pd(fix2,tx);
3047             fiy2             = _mm_add_pd(fiy2,ty);
3048             fiz2             = _mm_add_pd(fiz2,tz);
3049
3050             fjx2             = _mm_add_pd(fjx2,tx);
3051             fjy2             = _mm_add_pd(fjy2,ty);
3052             fjz2             = _mm_add_pd(fjz2,tz);
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 573 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*573);
3078 }