src/gmxlib/nonbonded/nb_kernel_sse4_1_single/nb_kernel_ElecRFCut_VdwLJSw_GeomW3W3_sse4_1_single.c

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