src/gmxlib/nonbonded/nb_kernel_avx_128_fma_single/nb_kernel_ElecEwSw_VdwNone_GeomW3W3_avx_128_fma_single.c

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