src/gmxlib/nonbonded/nb_kernel_sse4_1_double/nb_kernel_ElecEwSh_VdwLJSh_GeomW4W4_sse4_1_double.c

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