src/gromacs/gmxlib/nonbonded/nb_kernel_sse2_single/nb_kernel_ElecEwSw_VdwNone_GeomW4W4_sse2_single.c

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