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

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