src/gromacs/gmxlib/nonbonded/nb_kernel_sse4_1_double/nb_kernel_ElecEwSh_VdwNone_GeomW4W4_sse4_1_double.cpp

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