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

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