src/gromacs/gmxlib/nonbonded/nb_kernel_avx_128_fma_single/nb_kernel_ElecRF_VdwLJ_GeomW4W4_avx_128_fma_single.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*
  36  * Note: this file was generated by the GROMACS avx_128_fma_single kernel generator.
  37  */
  38 #include "config.h"
  39
  40 #include <math.h>
  41
  42 #include "../nb_kernel.h"
  43 #include "gromacs/legacyheaders/types/simple.h"
  44 #include "gromacs/math/vec.h"
  45 #include "gromacs/legacyheaders/nrnb.h"
  46
  47 #include "gromacs/simd/math_x86_avx_128_fma_single.h"
  48 #include "kernelutil_x86_avx_128_fma_single.h"
  49
  50 /*
  51  * Gromacs nonbonded kernel:   nb_kernel_ElecRF_VdwLJ_GeomW4W4_VF_avx_128_fma_single
  52  * Electrostatics interaction: ReactionField
  53  * VdW interaction:            LennardJones
  54  * Geometry:                   Water4-Water4
  55  * Calculate force/pot:        PotentialAndForce
  56  */
  57 void
  58 nb_kernel_ElecRF_VdwLJ_GeomW4W4_VF_avx_128_fma_single
  59                     (t_nblist                    * gmx_restrict       nlist,
  60                      rvec                        * gmx_restrict          xx,
  61                      rvec                        * gmx_restrict          ff,
  62                      t_forcerec                  * gmx_restrict          fr,
  63                      t_mdatoms                   * gmx_restrict     mdatoms,
  64                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
  65                      t_nrnb                      * gmx_restrict        nrnb)
  66 {
  67     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
  68      * just 0 for non-waters.
  69      * Suffixes A,B,C,D refer to j loop unrolling done with AVX_128, e.g. for the four different
  70      * jnr indices corresponding to data put in the four positions in the SIMD register.
  71      */
  72     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
  73     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
  74     int              jnrA,jnrB,jnrC,jnrD;
  75     int              jnrlistA,jnrlistB,jnrlistC,jnrlistD;
  76     int              j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
  77     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
  78     real             rcutoff_scalar;
  79     real             *shiftvec,*fshift,*x,*f;
  80     real             *fjptrA,*fjptrB,*fjptrC,*fjptrD;
  81     real             scratch[4*DIM];
  82     __m128           fscal,rcutoff,rcutoff2,jidxall;
  83     int              vdwioffset0;
  84     __m128           ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
  85     int              vdwioffset1;
  86     __m128           ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
  87     int              vdwioffset2;
  88     __m128           ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
  89     int              vdwioffset3;
  90     __m128           ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
  91     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
  92     __m128           jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
  93     int              vdwjidx1A,vdwjidx1B,vdwjidx1C,vdwjidx1D;
  94     __m128           jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
  95     int              vdwjidx2A,vdwjidx2B,vdwjidx2C,vdwjidx2D;
  96     __m128           jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
  97     int              vdwjidx3A,vdwjidx3B,vdwjidx3C,vdwjidx3D;
  98     __m128           jx3,jy3,jz3,fjx3,fjy3,fjz3,jq3,isaj3;
  99     __m128           dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
 100     __m128           dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11;
 101     __m128           dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12;
 102     __m128           dx13,dy13,dz13,rsq13,rinv13,rinvsq13,r13,qq13,c6_13,c12_13;
 103     __m128           dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21;
 104     __m128           dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22;
 105     __m128           dx23,dy23,dz23,rsq23,rinv23,rinvsq23,r23,qq23,c6_23,c12_23;
 106     __m128           dx31,dy31,dz31,rsq31,rinv31,rinvsq31,r31,qq31,c6_31,c12_31;
 107     __m128           dx32,dy32,dz32,rsq32,rinv32,rinvsq32,r32,qq32,c6_32,c12_32;
 108     __m128           dx33,dy33,dz33,rsq33,rinv33,rinvsq33,r33,qq33,c6_33,c12_33;
 109     __m128           velec,felec,velecsum,facel,crf,krf,krf2;
 110     real             *charge;
 111     int              nvdwtype;
 112     __m128           rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
 113     int              *vdwtype;
 114     real             *vdwparam;
 115     __m128           one_sixth   = _mm_set1_ps(1.0/6.0);
 116     __m128           one_twelfth = _mm_set1_ps(1.0/12.0);
 117     __m128           dummy_mask,cutoff_mask;
 118     __m128           signbit = _mm_castsi128_ps( _mm_set1_epi32(0x80000000) );
 119     __m128           one     = _mm_set1_ps(1.0);
 120     __m128           two     = _mm_set1_ps(2.0);
 121     x                = xx[0];
 122     f                = ff[0];
 123
 124     nri              = nlist->nri;
 125     iinr             = nlist->iinr;
 126     jindex           = nlist->jindex;
 127     jjnr             = nlist->jjnr;
 128     shiftidx         = nlist->shift;
 129     gid              = nlist->gid;
 130     shiftvec         = fr->shift_vec[0];
 131     fshift           = fr->fshift[0];
 132     facel            = _mm_set1_ps(fr->epsfac);
 133     charge           = mdatoms->chargeA;
 134     krf              = _mm_set1_ps(fr->ic->k_rf);
 135     krf2             = _mm_set1_ps(fr->ic->k_rf*2.0);
 136     crf              = _mm_set1_ps(fr->ic->c_rf);
 137     nvdwtype         = fr->ntype;
 138     vdwparam         = fr->nbfp;
 139     vdwtype          = mdatoms->typeA;
 140
 141     /* Setup water-specific parameters */
 142     inr              = nlist->iinr[0];
 143     iq1              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+1]));
 144     iq2              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+2]));
 145     iq3              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+3]));
 146     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
 147
 148     jq1              = _mm_set1_ps(charge[inr+1]);
 149     jq2              = _mm_set1_ps(charge[inr+2]);
 150     jq3              = _mm_set1_ps(charge[inr+3]);
 151     vdwjidx0A        = 2*vdwtype[inr+0];
 152     c6_00            = _mm_set1_ps(vdwparam[vdwioffset0+vdwjidx0A]);
 153     c12_00           = _mm_set1_ps(vdwparam[vdwioffset0+vdwjidx0A+1]);
 154     qq11             = _mm_mul_ps(iq1,jq1);
 155     qq12             = _mm_mul_ps(iq1,jq2);
 156     qq13             = _mm_mul_ps(iq1,jq3);
 157     qq21             = _mm_mul_ps(iq2,jq1);
 158     qq22             = _mm_mul_ps(iq2,jq2);
 159     qq23             = _mm_mul_ps(iq2,jq3);
 160     qq31             = _mm_mul_ps(iq3,jq1);
 161     qq32             = _mm_mul_ps(iq3,jq2);
 162     qq33             = _mm_mul_ps(iq3,jq3);
 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_4rvec_broadcast_ps(shiftvec+i_shift_offset,x+i_coord_offset,
 195                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
 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         fix3             = _mm_setzero_ps();
 207         fiy3             = _mm_setzero_ps();
 208         fiz3             = _mm_setzero_ps();
 209
 210         /* Reset potential sums */
 211         velecsum         = _mm_setzero_ps();
 212         vvdwsum          = _mm_setzero_ps();
 213
 214         /* Start inner kernel loop */
 215         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
 216         {
 217
 218             /* Get j neighbor index, and coordinate index */
 219             jnrA             = jjnr[jidx];
 220             jnrB             = jjnr[jidx+1];
 221             jnrC             = jjnr[jidx+2];
 222             jnrD             = jjnr[jidx+3];
 223             j_coord_offsetA  = DIM*jnrA;
 224             j_coord_offsetB  = DIM*jnrB;
 225             j_coord_offsetC  = DIM*jnrC;
 226             j_coord_offsetD  = DIM*jnrD;
 227
 228             /* load j atom coordinates */
 229             gmx_mm_load_4rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
 230                                               x+j_coord_offsetC,x+j_coord_offsetD,
 231                                               &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,
 232                                               &jy2,&jz2,&jx3,&jy3,&jz3);
 233
 234             /* Calculate displacement vector */
 235             dx00             = _mm_sub_ps(ix0,jx0);
 236             dy00             = _mm_sub_ps(iy0,jy0);
 237             dz00             = _mm_sub_ps(iz0,jz0);
 238             dx11             = _mm_sub_ps(ix1,jx1);
 239             dy11             = _mm_sub_ps(iy1,jy1);
 240             dz11             = _mm_sub_ps(iz1,jz1);
 241             dx12             = _mm_sub_ps(ix1,jx2);
 242             dy12             = _mm_sub_ps(iy1,jy2);
 243             dz12             = _mm_sub_ps(iz1,jz2);
 244             dx13             = _mm_sub_ps(ix1,jx3);
 245             dy13             = _mm_sub_ps(iy1,jy3);
 246             dz13             = _mm_sub_ps(iz1,jz3);
 247             dx21             = _mm_sub_ps(ix2,jx1);
 248             dy21             = _mm_sub_ps(iy2,jy1);
 249             dz21             = _mm_sub_ps(iz2,jz1);
 250             dx22             = _mm_sub_ps(ix2,jx2);
 251             dy22             = _mm_sub_ps(iy2,jy2);
 252             dz22             = _mm_sub_ps(iz2,jz2);
 253             dx23             = _mm_sub_ps(ix2,jx3);
 254             dy23             = _mm_sub_ps(iy2,jy3);
 255             dz23             = _mm_sub_ps(iz2,jz3);
 256             dx31             = _mm_sub_ps(ix3,jx1);
 257             dy31             = _mm_sub_ps(iy3,jy1);
 258             dz31             = _mm_sub_ps(iz3,jz1);
 259             dx32             = _mm_sub_ps(ix3,jx2);
 260             dy32             = _mm_sub_ps(iy3,jy2);
 261             dz32             = _mm_sub_ps(iz3,jz2);
 262             dx33             = _mm_sub_ps(ix3,jx3);
 263             dy33             = _mm_sub_ps(iy3,jy3);
 264             dz33             = _mm_sub_ps(iz3,jz3);
 265
 266             /* Calculate squared distance and things based on it */
 267             rsq00            = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
 268             rsq11            = gmx_mm_calc_rsq_ps(dx11,dy11,dz11);
 269             rsq12            = gmx_mm_calc_rsq_ps(dx12,dy12,dz12);
 270             rsq13            = gmx_mm_calc_rsq_ps(dx13,dy13,dz13);
 271             rsq21            = gmx_mm_calc_rsq_ps(dx21,dy21,dz21);
 272             rsq22            = gmx_mm_calc_rsq_ps(dx22,dy22,dz22);
 273             rsq23            = gmx_mm_calc_rsq_ps(dx23,dy23,dz23);
 274             rsq31            = gmx_mm_calc_rsq_ps(dx31,dy31,dz31);
 275             rsq32            = gmx_mm_calc_rsq_ps(dx32,dy32,dz32);
 276             rsq33            = gmx_mm_calc_rsq_ps(dx33,dy33,dz33);
 277
 278             rinv11           = gmx_mm_invsqrt_ps(rsq11);
 279             rinv12           = gmx_mm_invsqrt_ps(rsq12);
 280             rinv13           = gmx_mm_invsqrt_ps(rsq13);
 281             rinv21           = gmx_mm_invsqrt_ps(rsq21);
 282             rinv22           = gmx_mm_invsqrt_ps(rsq22);
 283             rinv23           = gmx_mm_invsqrt_ps(rsq23);
 284             rinv31           = gmx_mm_invsqrt_ps(rsq31);
 285             rinv32           = gmx_mm_invsqrt_ps(rsq32);
 286             rinv33           = gmx_mm_invsqrt_ps(rsq33);
 287
 288             rinvsq00         = gmx_mm_inv_ps(rsq00);
 289             rinvsq11         = _mm_mul_ps(rinv11,rinv11);
 290             rinvsq12         = _mm_mul_ps(rinv12,rinv12);
 291             rinvsq13         = _mm_mul_ps(rinv13,rinv13);
 292             rinvsq21         = _mm_mul_ps(rinv21,rinv21);
 293             rinvsq22         = _mm_mul_ps(rinv22,rinv22);
 294             rinvsq23         = _mm_mul_ps(rinv23,rinv23);
 295             rinvsq31         = _mm_mul_ps(rinv31,rinv31);
 296             rinvsq32         = _mm_mul_ps(rinv32,rinv32);
 297             rinvsq33         = _mm_mul_ps(rinv33,rinv33);
 298
 299             fjx0             = _mm_setzero_ps();
 300             fjy0             = _mm_setzero_ps();
 301             fjz0             = _mm_setzero_ps();
 302             fjx1             = _mm_setzero_ps();
 303             fjy1             = _mm_setzero_ps();
 304             fjz1             = _mm_setzero_ps();
 305             fjx2             = _mm_setzero_ps();
 306             fjy2             = _mm_setzero_ps();
 307             fjz2             = _mm_setzero_ps();
 308             fjx3             = _mm_setzero_ps();
 309             fjy3             = _mm_setzero_ps();
 310             fjz3             = _mm_setzero_ps();
 311
 312             /**************************
 313              * CALCULATE INTERACTIONS *
 314              **************************/
 315
 316             /* LENNARD-JONES DISPERSION/REPULSION */
 317
 318             rinvsix          = _mm_mul_ps(_mm_mul_ps(rinvsq00,rinvsq00),rinvsq00);
 319             vvdw6            = _mm_mul_ps(c6_00,rinvsix);
 320             vvdw12           = _mm_mul_ps(c12_00,_mm_mul_ps(rinvsix,rinvsix));
 321             vvdw             = _mm_msub_ps(vvdw12,one_twelfth,_mm_mul_ps(vvdw6,one_sixth));
 322             fvdw             = _mm_mul_ps(_mm_sub_ps(vvdw12,vvdw6),rinvsq00);
 323
 324             /* Update potential sum for this i atom from the interaction with this j atom. */
 325             vvdwsum          = _mm_add_ps(vvdwsum,vvdw);
 326
 327             fscal            = fvdw;
 328
 329              /* Update vectorial force */
 330             fix0             = _mm_macc_ps(dx00,fscal,fix0);
 331             fiy0             = _mm_macc_ps(dy00,fscal,fiy0);
 332             fiz0             = _mm_macc_ps(dz00,fscal,fiz0);
 333
 334             fjx0             = _mm_macc_ps(dx00,fscal,fjx0);
 335             fjy0             = _mm_macc_ps(dy00,fscal,fjy0);
 336             fjz0             = _mm_macc_ps(dz00,fscal,fjz0);
 337
 338             /**************************
 339              * CALCULATE INTERACTIONS *
 340              **************************/
 341
 342             /* REACTION-FIELD ELECTROSTATICS */
 343             velec            = _mm_mul_ps(qq11,_mm_sub_ps(_mm_macc_ps(krf,rsq11,rinv11),crf));
 344             felec            = _mm_mul_ps(qq11,_mm_msub_ps(rinv11,rinvsq11,krf2));
 345
 346             /* Update potential sum for this i atom from the interaction with this j atom. */
 347             velecsum         = _mm_add_ps(velecsum,velec);
 348
 349             fscal            = felec;
 350
 351              /* Update vectorial force */
 352             fix1             = _mm_macc_ps(dx11,fscal,fix1);
 353             fiy1             = _mm_macc_ps(dy11,fscal,fiy1);
 354             fiz1             = _mm_macc_ps(dz11,fscal,fiz1);
 355
 356             fjx1             = _mm_macc_ps(dx11,fscal,fjx1);
 357             fjy1             = _mm_macc_ps(dy11,fscal,fjy1);
 358             fjz1             = _mm_macc_ps(dz11,fscal,fjz1);
 359
 360             /**************************
 361              * CALCULATE INTERACTIONS *
 362              **************************/
 363
 364             /* REACTION-FIELD ELECTROSTATICS */
 365             velec            = _mm_mul_ps(qq12,_mm_sub_ps(_mm_macc_ps(krf,rsq12,rinv12),crf));
 366             felec            = _mm_mul_ps(qq12,_mm_msub_ps(rinv12,rinvsq12,krf2));
 367
 368             /* Update potential sum for this i atom from the interaction with this j atom. */
 369             velecsum         = _mm_add_ps(velecsum,velec);
 370
 371             fscal            = felec;
 372
 373              /* Update vectorial force */
 374             fix1             = _mm_macc_ps(dx12,fscal,fix1);
 375             fiy1             = _mm_macc_ps(dy12,fscal,fiy1);
 376             fiz1             = _mm_macc_ps(dz12,fscal,fiz1);
 377
 378             fjx2             = _mm_macc_ps(dx12,fscal,fjx2);
 379             fjy2             = _mm_macc_ps(dy12,fscal,fjy2);
 380             fjz2             = _mm_macc_ps(dz12,fscal,fjz2);
 381
 382             /**************************
 383              * CALCULATE INTERACTIONS *
 384              **************************/
 385
 386             /* REACTION-FIELD ELECTROSTATICS */
 387             velec            = _mm_mul_ps(qq13,_mm_sub_ps(_mm_macc_ps(krf,rsq13,rinv13),crf));
 388             felec            = _mm_mul_ps(qq13,_mm_msub_ps(rinv13,rinvsq13,krf2));
 389
 390             /* Update potential sum for this i atom from the interaction with this j atom. */
 391             velecsum         = _mm_add_ps(velecsum,velec);
 392
 393             fscal            = felec;
 394
 395              /* Update vectorial force */
 396             fix1             = _mm_macc_ps(dx13,fscal,fix1);
 397             fiy1             = _mm_macc_ps(dy13,fscal,fiy1);
 398             fiz1             = _mm_macc_ps(dz13,fscal,fiz1);
 399
 400             fjx3             = _mm_macc_ps(dx13,fscal,fjx3);
 401             fjy3             = _mm_macc_ps(dy13,fscal,fjy3);
 402             fjz3             = _mm_macc_ps(dz13,fscal,fjz3);
 403
 404             /**************************
 405              * CALCULATE INTERACTIONS *
 406              **************************/
 407
 408             /* REACTION-FIELD ELECTROSTATICS */
 409             velec            = _mm_mul_ps(qq21,_mm_sub_ps(_mm_macc_ps(krf,rsq21,rinv21),crf));
 410             felec            = _mm_mul_ps(qq21,_mm_msub_ps(rinv21,rinvsq21,krf2));
 411
 412             /* Update potential sum for this i atom from the interaction with this j atom. */
 413             velecsum         = _mm_add_ps(velecsum,velec);
 414
 415             fscal            = felec;
 416
 417              /* Update vectorial force */
 418             fix2             = _mm_macc_ps(dx21,fscal,fix2);
 419             fiy2             = _mm_macc_ps(dy21,fscal,fiy2);
 420             fiz2             = _mm_macc_ps(dz21,fscal,fiz2);
 421
 422             fjx1             = _mm_macc_ps(dx21,fscal,fjx1);
 423             fjy1             = _mm_macc_ps(dy21,fscal,fjy1);
 424             fjz1             = _mm_macc_ps(dz21,fscal,fjz1);
 425
 426             /**************************
 427              * CALCULATE INTERACTIONS *
 428              **************************/
 429
 430             /* REACTION-FIELD ELECTROSTATICS */
 431             velec            = _mm_mul_ps(qq22,_mm_sub_ps(_mm_macc_ps(krf,rsq22,rinv22),crf));
 432             felec            = _mm_mul_ps(qq22,_mm_msub_ps(rinv22,rinvsq22,krf2));
 433
 434             /* Update potential sum for this i atom from the interaction with this j atom. */
 435             velecsum         = _mm_add_ps(velecsum,velec);
 436
 437             fscal            = felec;
 438
 439              /* Update vectorial force */
 440             fix2             = _mm_macc_ps(dx22,fscal,fix2);
 441             fiy2             = _mm_macc_ps(dy22,fscal,fiy2);
 442             fiz2             = _mm_macc_ps(dz22,fscal,fiz2);
 443
 444             fjx2             = _mm_macc_ps(dx22,fscal,fjx2);
 445             fjy2             = _mm_macc_ps(dy22,fscal,fjy2);
 446             fjz2             = _mm_macc_ps(dz22,fscal,fjz2);
 447
 448             /**************************
 449              * CALCULATE INTERACTIONS *
 450              **************************/
 451
 452             /* REACTION-FIELD ELECTROSTATICS */
 453             velec            = _mm_mul_ps(qq23,_mm_sub_ps(_mm_macc_ps(krf,rsq23,rinv23),crf));
 454             felec            = _mm_mul_ps(qq23,_mm_msub_ps(rinv23,rinvsq23,krf2));
 455
 456             /* Update potential sum for this i atom from the interaction with this j atom. */
 457             velecsum         = _mm_add_ps(velecsum,velec);
 458
 459             fscal            = felec;
 460
 461              /* Update vectorial force */
 462             fix2             = _mm_macc_ps(dx23,fscal,fix2);
 463             fiy2             = _mm_macc_ps(dy23,fscal,fiy2);
 464             fiz2             = _mm_macc_ps(dz23,fscal,fiz2);
 465
 466             fjx3             = _mm_macc_ps(dx23,fscal,fjx3);
 467             fjy3             = _mm_macc_ps(dy23,fscal,fjy3);
 468             fjz3             = _mm_macc_ps(dz23,fscal,fjz3);
 469
 470             /**************************
 471              * CALCULATE INTERACTIONS *
 472              **************************/
 473
 474             /* REACTION-FIELD ELECTROSTATICS */
 475             velec            = _mm_mul_ps(qq31,_mm_sub_ps(_mm_macc_ps(krf,rsq31,rinv31),crf));
 476             felec            = _mm_mul_ps(qq31,_mm_msub_ps(rinv31,rinvsq31,krf2));
 477
 478             /* Update potential sum for this i atom from the interaction with this j atom. */
 479             velecsum         = _mm_add_ps(velecsum,velec);
 480
 481             fscal            = felec;
 482
 483              /* Update vectorial force */
 484             fix3             = _mm_macc_ps(dx31,fscal,fix3);
 485             fiy3             = _mm_macc_ps(dy31,fscal,fiy3);
 486             fiz3             = _mm_macc_ps(dz31,fscal,fiz3);
 487
 488             fjx1             = _mm_macc_ps(dx31,fscal,fjx1);
 489             fjy1             = _mm_macc_ps(dy31,fscal,fjy1);
 490             fjz1             = _mm_macc_ps(dz31,fscal,fjz1);
 491
 492             /**************************
 493              * CALCULATE INTERACTIONS *
 494              **************************/
 495
 496             /* REACTION-FIELD ELECTROSTATICS */
 497             velec            = _mm_mul_ps(qq32,_mm_sub_ps(_mm_macc_ps(krf,rsq32,rinv32),crf));
 498             felec            = _mm_mul_ps(qq32,_mm_msub_ps(rinv32,rinvsq32,krf2));
 499
 500             /* Update potential sum for this i atom from the interaction with this j atom. */
 501             velecsum         = _mm_add_ps(velecsum,velec);
 502
 503             fscal            = felec;
 504
 505              /* Update vectorial force */
 506             fix3             = _mm_macc_ps(dx32,fscal,fix3);
 507             fiy3             = _mm_macc_ps(dy32,fscal,fiy3);
 508             fiz3             = _mm_macc_ps(dz32,fscal,fiz3);
 509
 510             fjx2             = _mm_macc_ps(dx32,fscal,fjx2);
 511             fjy2             = _mm_macc_ps(dy32,fscal,fjy2);
 512             fjz2             = _mm_macc_ps(dz32,fscal,fjz2);
 513
 514             /**************************
 515              * CALCULATE INTERACTIONS *
 516              **************************/
 517
 518             /* REACTION-FIELD ELECTROSTATICS */
 519             velec            = _mm_mul_ps(qq33,_mm_sub_ps(_mm_macc_ps(krf,rsq33,rinv33),crf));
 520             felec            = _mm_mul_ps(qq33,_mm_msub_ps(rinv33,rinvsq33,krf2));
 521
 522             /* Update potential sum for this i atom from the interaction with this j atom. */
 523             velecsum         = _mm_add_ps(velecsum,velec);
 524
 525             fscal            = felec;
 526
 527              /* Update vectorial force */
 528             fix3             = _mm_macc_ps(dx33,fscal,fix3);
 529             fiy3             = _mm_macc_ps(dy33,fscal,fiy3);
 530             fiz3             = _mm_macc_ps(dz33,fscal,fiz3);
 531
 532             fjx3             = _mm_macc_ps(dx33,fscal,fjx3);
 533             fjy3             = _mm_macc_ps(dy33,fscal,fjy3);
 534             fjz3             = _mm_macc_ps(dz33,fscal,fjz3);
 535
 536             fjptrA             = f+j_coord_offsetA;
 537             fjptrB             = f+j_coord_offsetB;
 538             fjptrC             = f+j_coord_offsetC;
 539             fjptrD             = f+j_coord_offsetD;
 540
 541             gmx_mm_decrement_4rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,
 542                                                    fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,
 543                                                    fjx2,fjy2,fjz2,fjx3,fjy3,fjz3);
 544
 545             /* Inner loop uses 353 flops */
 546         }
 547
 548         if(jidx<j_index_end)
 549         {
 550
 551             /* Get j neighbor index, and coordinate index */
 552             jnrlistA         = jjnr[jidx];
 553             jnrlistB         = jjnr[jidx+1];
 554             jnrlistC         = jjnr[jidx+2];
 555             jnrlistD         = jjnr[jidx+3];
 556             /* Sign of each element will be negative for non-real atoms.
 557              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
 558              * so use it as val = _mm_andnot_ps(mask,val) to clear dummy entries.
 559              */
 560             dummy_mask = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
 561             jnrA       = (jnrlistA>=0) ? jnrlistA : 0;
 562             jnrB       = (jnrlistB>=0) ? jnrlistB : 0;
 563             jnrC       = (jnrlistC>=0) ? jnrlistC : 0;
 564             jnrD       = (jnrlistD>=0) ? jnrlistD : 0;
 565             j_coord_offsetA  = DIM*jnrA;
 566             j_coord_offsetB  = DIM*jnrB;
 567             j_coord_offsetC  = DIM*jnrC;
 568             j_coord_offsetD  = DIM*jnrD;
 569
 570             /* load j atom coordinates */
 571             gmx_mm_load_4rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
 572                                               x+j_coord_offsetC,x+j_coord_offsetD,
 573                                               &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,
 574                                               &jy2,&jz2,&jx3,&jy3,&jz3);
 575
 576             /* Calculate displacement vector */
 577             dx00             = _mm_sub_ps(ix0,jx0);
 578             dy00             = _mm_sub_ps(iy0,jy0);
 579             dz00             = _mm_sub_ps(iz0,jz0);
 580             dx11             = _mm_sub_ps(ix1,jx1);
 581             dy11             = _mm_sub_ps(iy1,jy1);
 582             dz11             = _mm_sub_ps(iz1,jz1);
 583             dx12             = _mm_sub_ps(ix1,jx2);
 584             dy12             = _mm_sub_ps(iy1,jy2);
 585             dz12             = _mm_sub_ps(iz1,jz2);
 586             dx13             = _mm_sub_ps(ix1,jx3);
 587             dy13             = _mm_sub_ps(iy1,jy3);
 588             dz13             = _mm_sub_ps(iz1,jz3);
 589             dx21             = _mm_sub_ps(ix2,jx1);
 590             dy21             = _mm_sub_ps(iy2,jy1);
 591             dz21             = _mm_sub_ps(iz2,jz1);
 592             dx22             = _mm_sub_ps(ix2,jx2);
 593             dy22             = _mm_sub_ps(iy2,jy2);
 594             dz22             = _mm_sub_ps(iz2,jz2);
 595             dx23             = _mm_sub_ps(ix2,jx3);
 596             dy23             = _mm_sub_ps(iy2,jy3);
 597             dz23             = _mm_sub_ps(iz2,jz3);
 598             dx31             = _mm_sub_ps(ix3,jx1);
 599             dy31             = _mm_sub_ps(iy3,jy1);
 600             dz31             = _mm_sub_ps(iz3,jz1);
 601             dx32             = _mm_sub_ps(ix3,jx2);
 602             dy32             = _mm_sub_ps(iy3,jy2);
 603             dz32             = _mm_sub_ps(iz3,jz2);
 604             dx33             = _mm_sub_ps(ix3,jx3);
 605             dy33             = _mm_sub_ps(iy3,jy3);
 606             dz33             = _mm_sub_ps(iz3,jz3);
 607
 608             /* Calculate squared distance and things based on it */
 609             rsq00            = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
 610             rsq11            = gmx_mm_calc_rsq_ps(dx11,dy11,dz11);
 611             rsq12            = gmx_mm_calc_rsq_ps(dx12,dy12,dz12);
 612             rsq13            = gmx_mm_calc_rsq_ps(dx13,dy13,dz13);
 613             rsq21            = gmx_mm_calc_rsq_ps(dx21,dy21,dz21);
 614             rsq22            = gmx_mm_calc_rsq_ps(dx22,dy22,dz22);
 615             rsq23            = gmx_mm_calc_rsq_ps(dx23,dy23,dz23);
 616             rsq31            = gmx_mm_calc_rsq_ps(dx31,dy31,dz31);
 617             rsq32            = gmx_mm_calc_rsq_ps(dx32,dy32,dz32);
 618             rsq33            = gmx_mm_calc_rsq_ps(dx33,dy33,dz33);
 619
 620             rinv11           = gmx_mm_invsqrt_ps(rsq11);
 621             rinv12           = gmx_mm_invsqrt_ps(rsq12);
 622             rinv13           = gmx_mm_invsqrt_ps(rsq13);
 623             rinv21           = gmx_mm_invsqrt_ps(rsq21);
 624             rinv22           = gmx_mm_invsqrt_ps(rsq22);
 625             rinv23           = gmx_mm_invsqrt_ps(rsq23);
 626             rinv31           = gmx_mm_invsqrt_ps(rsq31);
 627             rinv32           = gmx_mm_invsqrt_ps(rsq32);
 628             rinv33           = gmx_mm_invsqrt_ps(rsq33);
 629
 630             rinvsq00         = gmx_mm_inv_ps(rsq00);
 631             rinvsq11         = _mm_mul_ps(rinv11,rinv11);
 632             rinvsq12         = _mm_mul_ps(rinv12,rinv12);
 633             rinvsq13         = _mm_mul_ps(rinv13,rinv13);
 634             rinvsq21         = _mm_mul_ps(rinv21,rinv21);
 635             rinvsq22         = _mm_mul_ps(rinv22,rinv22);
 636             rinvsq23         = _mm_mul_ps(rinv23,rinv23);
 637             rinvsq31         = _mm_mul_ps(rinv31,rinv31);
 638             rinvsq32         = _mm_mul_ps(rinv32,rinv32);
 639             rinvsq33         = _mm_mul_ps(rinv33,rinv33);
 640
 641             fjx0             = _mm_setzero_ps();
 642             fjy0             = _mm_setzero_ps();
 643             fjz0             = _mm_setzero_ps();
 644             fjx1             = _mm_setzero_ps();
 645             fjy1             = _mm_setzero_ps();
 646             fjz1             = _mm_setzero_ps();
 647             fjx2             = _mm_setzero_ps();
 648             fjy2             = _mm_setzero_ps();
 649             fjz2             = _mm_setzero_ps();
 650             fjx3             = _mm_setzero_ps();
 651             fjy3             = _mm_setzero_ps();
 652             fjz3             = _mm_setzero_ps();
 653
 654             /**************************
 655              * CALCULATE INTERACTIONS *
 656              **************************/
 657
 658             /* LENNARD-JONES DISPERSION/REPULSION */
 659
 660             rinvsix          = _mm_mul_ps(_mm_mul_ps(rinvsq00,rinvsq00),rinvsq00);
 661             vvdw6            = _mm_mul_ps(c6_00,rinvsix);
 662             vvdw12           = _mm_mul_ps(c12_00,_mm_mul_ps(rinvsix,rinvsix));
 663             vvdw             = _mm_msub_ps(vvdw12,one_twelfth,_mm_mul_ps(vvdw6,one_sixth));
 664             fvdw             = _mm_mul_ps(_mm_sub_ps(vvdw12,vvdw6),rinvsq00);
 665
 666             /* Update potential sum for this i atom from the interaction with this j atom. */
 667             vvdw             = _mm_andnot_ps(dummy_mask,vvdw);
 668             vvdwsum          = _mm_add_ps(vvdwsum,vvdw);
 669
 670             fscal            = fvdw;
 671
 672             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 673
 674              /* Update vectorial force */
 675             fix0             = _mm_macc_ps(dx00,fscal,fix0);
 676             fiy0             = _mm_macc_ps(dy00,fscal,fiy0);
 677             fiz0             = _mm_macc_ps(dz00,fscal,fiz0);
 678
 679             fjx0             = _mm_macc_ps(dx00,fscal,fjx0);
 680             fjy0             = _mm_macc_ps(dy00,fscal,fjy0);
 681             fjz0             = _mm_macc_ps(dz00,fscal,fjz0);
 682
 683             /**************************
 684              * CALCULATE INTERACTIONS *
 685              **************************/
 686
 687             /* REACTION-FIELD ELECTROSTATICS */
 688             velec            = _mm_mul_ps(qq11,_mm_sub_ps(_mm_macc_ps(krf,rsq11,rinv11),crf));
 689             felec            = _mm_mul_ps(qq11,_mm_msub_ps(rinv11,rinvsq11,krf2));
 690
 691             /* Update potential sum for this i atom from the interaction with this j atom. */
 692             velec            = _mm_andnot_ps(dummy_mask,velec);
 693             velecsum         = _mm_add_ps(velecsum,velec);
 694
 695             fscal            = felec;
 696
 697             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 698
 699              /* Update vectorial force */
 700             fix1             = _mm_macc_ps(dx11,fscal,fix1);
 701             fiy1             = _mm_macc_ps(dy11,fscal,fiy1);
 702             fiz1             = _mm_macc_ps(dz11,fscal,fiz1);
 703
 704             fjx1             = _mm_macc_ps(dx11,fscal,fjx1);
 705             fjy1             = _mm_macc_ps(dy11,fscal,fjy1);
 706             fjz1             = _mm_macc_ps(dz11,fscal,fjz1);
 707
 708             /**************************
 709              * CALCULATE INTERACTIONS *
 710              **************************/
 711
 712             /* REACTION-FIELD ELECTROSTATICS */
 713             velec            = _mm_mul_ps(qq12,_mm_sub_ps(_mm_macc_ps(krf,rsq12,rinv12),crf));
 714             felec            = _mm_mul_ps(qq12,_mm_msub_ps(rinv12,rinvsq12,krf2));
 715
 716             /* Update potential sum for this i atom from the interaction with this j atom. */
 717             velec            = _mm_andnot_ps(dummy_mask,velec);
 718             velecsum         = _mm_add_ps(velecsum,velec);
 719
 720             fscal            = felec;
 721
 722             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 723
 724              /* Update vectorial force */
 725             fix1             = _mm_macc_ps(dx12,fscal,fix1);
 726             fiy1             = _mm_macc_ps(dy12,fscal,fiy1);
 727             fiz1             = _mm_macc_ps(dz12,fscal,fiz1);
 728
 729             fjx2             = _mm_macc_ps(dx12,fscal,fjx2);
 730             fjy2             = _mm_macc_ps(dy12,fscal,fjy2);
 731             fjz2             = _mm_macc_ps(dz12,fscal,fjz2);
 732
 733             /**************************
 734              * CALCULATE INTERACTIONS *
 735              **************************/
 736
 737             /* REACTION-FIELD ELECTROSTATICS */
 738             velec            = _mm_mul_ps(qq13,_mm_sub_ps(_mm_macc_ps(krf,rsq13,rinv13),crf));
 739             felec            = _mm_mul_ps(qq13,_mm_msub_ps(rinv13,rinvsq13,krf2));
 740
 741             /* Update potential sum for this i atom from the interaction with this j atom. */
 742             velec            = _mm_andnot_ps(dummy_mask,velec);
 743             velecsum         = _mm_add_ps(velecsum,velec);
 744
 745             fscal            = felec;
 746
 747             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 748
 749              /* Update vectorial force */
 750             fix1             = _mm_macc_ps(dx13,fscal,fix1);
 751             fiy1             = _mm_macc_ps(dy13,fscal,fiy1);
 752             fiz1             = _mm_macc_ps(dz13,fscal,fiz1);
 753
 754             fjx3             = _mm_macc_ps(dx13,fscal,fjx3);
 755             fjy3             = _mm_macc_ps(dy13,fscal,fjy3);
 756             fjz3             = _mm_macc_ps(dz13,fscal,fjz3);
 757
 758             /**************************
 759              * CALCULATE INTERACTIONS *
 760              **************************/
 761
 762             /* REACTION-FIELD ELECTROSTATICS */
 763             velec            = _mm_mul_ps(qq21,_mm_sub_ps(_mm_macc_ps(krf,rsq21,rinv21),crf));
 764             felec            = _mm_mul_ps(qq21,_mm_msub_ps(rinv21,rinvsq21,krf2));
 765
 766             /* Update potential sum for this i atom from the interaction with this j atom. */
 767             velec            = _mm_andnot_ps(dummy_mask,velec);
 768             velecsum         = _mm_add_ps(velecsum,velec);
 769
 770             fscal            = felec;
 771
 772             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 773
 774              /* Update vectorial force */
 775             fix2             = _mm_macc_ps(dx21,fscal,fix2);
 776             fiy2             = _mm_macc_ps(dy21,fscal,fiy2);
 777             fiz2             = _mm_macc_ps(dz21,fscal,fiz2);
 778
 779             fjx1             = _mm_macc_ps(dx21,fscal,fjx1);
 780             fjy1             = _mm_macc_ps(dy21,fscal,fjy1);
 781             fjz1             = _mm_macc_ps(dz21,fscal,fjz1);
 782
 783             /**************************
 784              * CALCULATE INTERACTIONS *
 785              **************************/
 786
 787             /* REACTION-FIELD ELECTROSTATICS */
 788             velec            = _mm_mul_ps(qq22,_mm_sub_ps(_mm_macc_ps(krf,rsq22,rinv22),crf));
 789             felec            = _mm_mul_ps(qq22,_mm_msub_ps(rinv22,rinvsq22,krf2));
 790
 791             /* Update potential sum for this i atom from the interaction with this j atom. */
 792             velec            = _mm_andnot_ps(dummy_mask,velec);
 793             velecsum         = _mm_add_ps(velecsum,velec);
 794
 795             fscal            = felec;
 796
 797             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 798
 799              /* Update vectorial force */
 800             fix2             = _mm_macc_ps(dx22,fscal,fix2);
 801             fiy2             = _mm_macc_ps(dy22,fscal,fiy2);
 802             fiz2             = _mm_macc_ps(dz22,fscal,fiz2);
 803
 804             fjx2             = _mm_macc_ps(dx22,fscal,fjx2);
 805             fjy2             = _mm_macc_ps(dy22,fscal,fjy2);
 806             fjz2             = _mm_macc_ps(dz22,fscal,fjz2);
 807
 808             /**************************
 809              * CALCULATE INTERACTIONS *
 810              **************************/
 811
 812             /* REACTION-FIELD ELECTROSTATICS */
 813             velec            = _mm_mul_ps(qq23,_mm_sub_ps(_mm_macc_ps(krf,rsq23,rinv23),crf));
 814             felec            = _mm_mul_ps(qq23,_mm_msub_ps(rinv23,rinvsq23,krf2));
 815
 816             /* Update potential sum for this i atom from the interaction with this j atom. */
 817             velec            = _mm_andnot_ps(dummy_mask,velec);
 818             velecsum         = _mm_add_ps(velecsum,velec);
 819
 820             fscal            = felec;
 821
 822             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 823
 824              /* Update vectorial force */
 825             fix2             = _mm_macc_ps(dx23,fscal,fix2);
 826             fiy2             = _mm_macc_ps(dy23,fscal,fiy2);
 827             fiz2             = _mm_macc_ps(dz23,fscal,fiz2);
 828
 829             fjx3             = _mm_macc_ps(dx23,fscal,fjx3);
 830             fjy3             = _mm_macc_ps(dy23,fscal,fjy3);
 831             fjz3             = _mm_macc_ps(dz23,fscal,fjz3);
 832
 833             /**************************
 834              * CALCULATE INTERACTIONS *
 835              **************************/
 836
 837             /* REACTION-FIELD ELECTROSTATICS */
 838             velec            = _mm_mul_ps(qq31,_mm_sub_ps(_mm_macc_ps(krf,rsq31,rinv31),crf));
 839             felec            = _mm_mul_ps(qq31,_mm_msub_ps(rinv31,rinvsq31,krf2));
 840
 841             /* Update potential sum for this i atom from the interaction with this j atom. */
 842             velec            = _mm_andnot_ps(dummy_mask,velec);
 843             velecsum         = _mm_add_ps(velecsum,velec);
 844
 845             fscal            = felec;
 846
 847             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 848
 849              /* Update vectorial force */
 850             fix3             = _mm_macc_ps(dx31,fscal,fix3);
 851             fiy3             = _mm_macc_ps(dy31,fscal,fiy3);
 852             fiz3             = _mm_macc_ps(dz31,fscal,fiz3);
 853
 854             fjx1             = _mm_macc_ps(dx31,fscal,fjx1);
 855             fjy1             = _mm_macc_ps(dy31,fscal,fjy1);
 856             fjz1             = _mm_macc_ps(dz31,fscal,fjz1);
 857
 858             /**************************
 859              * CALCULATE INTERACTIONS *
 860              **************************/
 861
 862             /* REACTION-FIELD ELECTROSTATICS */
 863             velec            = _mm_mul_ps(qq32,_mm_sub_ps(_mm_macc_ps(krf,rsq32,rinv32),crf));
 864             felec            = _mm_mul_ps(qq32,_mm_msub_ps(rinv32,rinvsq32,krf2));
 865
 866             /* Update potential sum for this i atom from the interaction with this j atom. */
 867             velec            = _mm_andnot_ps(dummy_mask,velec);
 868             velecsum         = _mm_add_ps(velecsum,velec);
 869
 870             fscal            = felec;
 871
 872             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 873
 874              /* Update vectorial force */
 875             fix3             = _mm_macc_ps(dx32,fscal,fix3);
 876             fiy3             = _mm_macc_ps(dy32,fscal,fiy3);
 877             fiz3             = _mm_macc_ps(dz32,fscal,fiz3);
 878
 879             fjx2             = _mm_macc_ps(dx32,fscal,fjx2);
 880             fjy2             = _mm_macc_ps(dy32,fscal,fjy2);
 881             fjz2             = _mm_macc_ps(dz32,fscal,fjz2);
 882
 883             /**************************
 884              * CALCULATE INTERACTIONS *
 885              **************************/
 886
 887             /* REACTION-FIELD ELECTROSTATICS */
 888             velec            = _mm_mul_ps(qq33,_mm_sub_ps(_mm_macc_ps(krf,rsq33,rinv33),crf));
 889             felec            = _mm_mul_ps(qq33,_mm_msub_ps(rinv33,rinvsq33,krf2));
 890
 891             /* Update potential sum for this i atom from the interaction with this j atom. */
 892             velec            = _mm_andnot_ps(dummy_mask,velec);
 893             velecsum         = _mm_add_ps(velecsum,velec);
 894
 895             fscal            = felec;
 896
 897             fscal            = _mm_andnot_ps(dummy_mask,fscal);
 898
 899              /* Update vectorial force */
 900             fix3             = _mm_macc_ps(dx33,fscal,fix3);
 901             fiy3             = _mm_macc_ps(dy33,fscal,fiy3);
 902             fiz3             = _mm_macc_ps(dz33,fscal,fiz3);
 903
 904             fjx3             = _mm_macc_ps(dx33,fscal,fjx3);
 905             fjy3             = _mm_macc_ps(dy33,fscal,fjy3);
 906             fjz3             = _mm_macc_ps(dz33,fscal,fjz3);
 907
 908             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
 909             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
 910             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
 911             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
 912
 913             gmx_mm_decrement_4rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,
 914                                                    fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,
 915                                                    fjx2,fjy2,fjz2,fjx3,fjy3,fjz3);
 916
 917             /* Inner loop uses 353 flops */
 918         }
 919
 920         /* End of innermost loop */
 921
 922         gmx_mm_update_iforce_4atom_swizzle_ps(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
 923                                               f+i_coord_offset,fshift+i_shift_offset);
 924
 925         ggid                        = gid[iidx];
 926         /* Update potential energies */
 927         gmx_mm_update_1pot_ps(velecsum,kernel_data->energygrp_elec+ggid);
 928         gmx_mm_update_1pot_ps(vvdwsum,kernel_data->energygrp_vdw+ggid);
 929
 930         /* Increment number of inner iterations */
 931         inneriter                  += j_index_end - j_index_start;
 932
 933         /* Outer loop uses 26 flops */
 934     }
 935
 936     /* Increment number of outer iterations */
 937     outeriter        += nri;
 938
 939     /* Update outer/inner flops */
 940
 941     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4W4_VF,outeriter*26 + inneriter*353);
 942 }
 943 /*
 944  * Gromacs nonbonded kernel:   nb_kernel_ElecRF_VdwLJ_GeomW4W4_F_avx_128_fma_single
 945  * Electrostatics interaction: ReactionField
 946  * VdW interaction:            LennardJones
 947  * Geometry:                   Water4-Water4
 948  * Calculate force/pot:        Force
 949  */
 950 void
 951 nb_kernel_ElecRF_VdwLJ_GeomW4W4_F_avx_128_fma_single
 952                     (t_nblist                    * gmx_restrict       nlist,
 953                      rvec                        * gmx_restrict          xx,
 954                      rvec                        * gmx_restrict          ff,
 955                      t_forcerec                  * gmx_restrict          fr,
 956                      t_mdatoms                   * gmx_restrict     mdatoms,
 957                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
 958                      t_nrnb                      * gmx_restrict        nrnb)
 959 {
 960     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
 961      * just 0 for non-waters.
 962      * Suffixes A,B,C,D refer to j loop unrolling done with AVX_128, e.g. for the four different
 963      * jnr indices corresponding to data put in the four positions in the SIMD register.
 964      */
 965     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
 966     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
 967     int              jnrA,jnrB,jnrC,jnrD;
 968     int              jnrlistA,jnrlistB,jnrlistC,jnrlistD;
 969     int              j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
 970     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
 971     real             rcutoff_scalar;
 972     real             *shiftvec,*fshift,*x,*f;
 973     real             *fjptrA,*fjptrB,*fjptrC,*fjptrD;
 974     real             scratch[4*DIM];
 975     __m128           fscal,rcutoff,rcutoff2,jidxall;
 976     int              vdwioffset0;
 977     __m128           ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
 978     int              vdwioffset1;
 979     __m128           ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
 980     int              vdwioffset2;
 981     __m128           ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
 982     int              vdwioffset3;
 983     __m128           ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
 984     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
 985     __m128           jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
 986     int              vdwjidx1A,vdwjidx1B,vdwjidx1C,vdwjidx1D;
 987     __m128           jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
 988     int              vdwjidx2A,vdwjidx2B,vdwjidx2C,vdwjidx2D;
 989     __m128           jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
 990     int              vdwjidx3A,vdwjidx3B,vdwjidx3C,vdwjidx3D;
 991     __m128           jx3,jy3,jz3,fjx3,fjy3,fjz3,jq3,isaj3;
 992     __m128           dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
 993     __m128           dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11;
 994     __m128           dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12;
 995     __m128           dx13,dy13,dz13,rsq13,rinv13,rinvsq13,r13,qq13,c6_13,c12_13;
 996     __m128           dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21;
 997     __m128           dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22;
 998     __m128           dx23,dy23,dz23,rsq23,rinv23,rinvsq23,r23,qq23,c6_23,c12_23;
 999     __m128           dx31,dy31,dz31,rsq31,rinv31,rinvsq31,r31,qq31,c6_31,c12_31;
1000     __m128           dx32,dy32,dz32,rsq32,rinv32,rinvsq32,r32,qq32,c6_32,c12_32;
1001     __m128           dx33,dy33,dz33,rsq33,rinv33,rinvsq33,r33,qq33,c6_33,c12_33;
1002     __m128           velec,felec,velecsum,facel,crf,krf,krf2;
1003     real             *charge;
1004     int              nvdwtype;
1005     __m128           rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
1006     int              *vdwtype;
1007     real             *vdwparam;
1008     __m128           one_sixth   = _mm_set1_ps(1.0/6.0);
1009     __m128           one_twelfth = _mm_set1_ps(1.0/12.0);
1010     __m128           dummy_mask,cutoff_mask;
1011     __m128           signbit = _mm_castsi128_ps( _mm_set1_epi32(0x80000000) );
1012     __m128           one     = _mm_set1_ps(1.0);
1013     __m128           two     = _mm_set1_ps(2.0);
1014     x                = xx[0];
1015     f                = ff[0];
1016
1017     nri              = nlist->nri;
1018     iinr             = nlist->iinr;
1019     jindex           = nlist->jindex;
1020     jjnr             = nlist->jjnr;
1021     shiftidx         = nlist->shift;
1022     gid              = nlist->gid;
1023     shiftvec         = fr->shift_vec[0];
1024     fshift           = fr->fshift[0];
1025     facel            = _mm_set1_ps(fr->epsfac);
1026     charge           = mdatoms->chargeA;
1027     krf              = _mm_set1_ps(fr->ic->k_rf);
1028     krf2             = _mm_set1_ps(fr->ic->k_rf*2.0);
1029     crf              = _mm_set1_ps(fr->ic->c_rf);
1030     nvdwtype         = fr->ntype;
1031     vdwparam         = fr->nbfp;
1032     vdwtype          = mdatoms->typeA;
1033
1034     /* Setup water-specific parameters */
1035     inr              = nlist->iinr[0];
1036     iq1              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+1]));
1037     iq2              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+2]));
1038     iq3              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+3]));
1039     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
1040
1041     jq1              = _mm_set1_ps(charge[inr+1]);
1042     jq2              = _mm_set1_ps(charge[inr+2]);
1043     jq3              = _mm_set1_ps(charge[inr+3]);
1044     vdwjidx0A        = 2*vdwtype[inr+0];
1045     c6_00            = _mm_set1_ps(vdwparam[vdwioffset0+vdwjidx0A]);
1046     c12_00           = _mm_set1_ps(vdwparam[vdwioffset0+vdwjidx0A+1]);
1047     qq11             = _mm_mul_ps(iq1,jq1);
1048     qq12             = _mm_mul_ps(iq1,jq2);
1049     qq13             = _mm_mul_ps(iq1,jq3);
1050     qq21             = _mm_mul_ps(iq2,jq1);
1051     qq22             = _mm_mul_ps(iq2,jq2);
1052     qq23             = _mm_mul_ps(iq2,jq3);
1053     qq31             = _mm_mul_ps(iq3,jq1);
1054     qq32             = _mm_mul_ps(iq3,jq2);
1055     qq33             = _mm_mul_ps(iq3,jq3);
1056
1057     /* Avoid stupid compiler warnings */
1058     jnrA = jnrB = jnrC = jnrD = 0;
1059     j_coord_offsetA = 0;
1060     j_coord_offsetB = 0;
1061     j_coord_offsetC = 0;
1062     j_coord_offsetD = 0;
1063
1064     outeriter        = 0;
1065     inneriter        = 0;
1066
1067     for(iidx=0;iidx<4*DIM;iidx++)
1068     {
1069         scratch[iidx] = 0.0;
1070     }
1071
1072     /* Start outer loop over neighborlists */
1073     for(iidx=0; iidx<nri; iidx++)
1074     {
1075         /* Load shift vector for this list */
1076         i_shift_offset   = DIM*shiftidx[iidx];
1077
1078         /* Load limits for loop over neighbors */
1079         j_index_start    = jindex[iidx];
1080         j_index_end      = jindex[iidx+1];
1081
1082         /* Get outer coordinate index */
1083         inr              = iinr[iidx];
1084         i_coord_offset   = DIM*inr;
1085
1086         /* Load i particle coords and add shift vector */
1087         gmx_mm_load_shift_and_4rvec_broadcast_ps(shiftvec+i_shift_offset,x+i_coord_offset,
1088                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
1089
1090         fix0             = _mm_setzero_ps();
1091         fiy0             = _mm_setzero_ps();
1092         fiz0             = _mm_setzero_ps();
1093         fix1             = _mm_setzero_ps();
1094         fiy1             = _mm_setzero_ps();
1095         fiz1             = _mm_setzero_ps();
1096         fix2             = _mm_setzero_ps();
1097         fiy2             = _mm_setzero_ps();
1098         fiz2             = _mm_setzero_ps();
1099         fix3             = _mm_setzero_ps();
1100         fiy3             = _mm_setzero_ps();
1101         fiz3             = _mm_setzero_ps();
1102
1103         /* Start inner kernel loop */
1104         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
1105         {
1106
1107             /* Get j neighbor index, and coordinate index */
1108             jnrA             = jjnr[jidx];
1109             jnrB             = jjnr[jidx+1];
1110             jnrC             = jjnr[jidx+2];
1111             jnrD             = jjnr[jidx+3];
1112             j_coord_offsetA  = DIM*jnrA;
1113             j_coord_offsetB  = DIM*jnrB;
1114             j_coord_offsetC  = DIM*jnrC;
1115             j_coord_offsetD  = DIM*jnrD;
1116
1117             /* load j atom coordinates */
1118             gmx_mm_load_4rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
1119                                               x+j_coord_offsetC,x+j_coord_offsetD,
1120                                               &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,
1121                                               &jy2,&jz2,&jx3,&jy3,&jz3);
1122
1123             /* Calculate displacement vector */
1124             dx00             = _mm_sub_ps(ix0,jx0);
1125             dy00             = _mm_sub_ps(iy0,jy0);
1126             dz00             = _mm_sub_ps(iz0,jz0);
1127             dx11             = _mm_sub_ps(ix1,jx1);
1128             dy11             = _mm_sub_ps(iy1,jy1);
1129             dz11             = _mm_sub_ps(iz1,jz1);
1130             dx12             = _mm_sub_ps(ix1,jx2);
1131             dy12             = _mm_sub_ps(iy1,jy2);
1132             dz12             = _mm_sub_ps(iz1,jz2);
1133             dx13             = _mm_sub_ps(ix1,jx3);
1134             dy13             = _mm_sub_ps(iy1,jy3);
1135             dz13             = _mm_sub_ps(iz1,jz3);
1136             dx21             = _mm_sub_ps(ix2,jx1);
1137             dy21             = _mm_sub_ps(iy2,jy1);
1138             dz21             = _mm_sub_ps(iz2,jz1);
1139             dx22             = _mm_sub_ps(ix2,jx2);
1140             dy22             = _mm_sub_ps(iy2,jy2);
1141             dz22             = _mm_sub_ps(iz2,jz2);
1142             dx23             = _mm_sub_ps(ix2,jx3);
1143             dy23             = _mm_sub_ps(iy2,jy3);
1144             dz23             = _mm_sub_ps(iz2,jz3);
1145             dx31             = _mm_sub_ps(ix3,jx1);
1146             dy31             = _mm_sub_ps(iy3,jy1);
1147             dz31             = _mm_sub_ps(iz3,jz1);
1148             dx32             = _mm_sub_ps(ix3,jx2);
1149             dy32             = _mm_sub_ps(iy3,jy2);
1150             dz32             = _mm_sub_ps(iz3,jz2);
1151             dx33             = _mm_sub_ps(ix3,jx3);
1152             dy33             = _mm_sub_ps(iy3,jy3);
1153             dz33             = _mm_sub_ps(iz3,jz3);
1154
1155             /* Calculate squared distance and things based on it */
1156             rsq00            = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
1157             rsq11            = gmx_mm_calc_rsq_ps(dx11,dy11,dz11);
1158             rsq12            = gmx_mm_calc_rsq_ps(dx12,dy12,dz12);
1159             rsq13            = gmx_mm_calc_rsq_ps(dx13,dy13,dz13);
1160             rsq21            = gmx_mm_calc_rsq_ps(dx21,dy21,dz21);
1161             rsq22            = gmx_mm_calc_rsq_ps(dx22,dy22,dz22);
1162             rsq23            = gmx_mm_calc_rsq_ps(dx23,dy23,dz23);
1163             rsq31            = gmx_mm_calc_rsq_ps(dx31,dy31,dz31);
1164             rsq32            = gmx_mm_calc_rsq_ps(dx32,dy32,dz32);
1165             rsq33            = gmx_mm_calc_rsq_ps(dx33,dy33,dz33);
1166
1167             rinv11           = gmx_mm_invsqrt_ps(rsq11);
1168             rinv12           = gmx_mm_invsqrt_ps(rsq12);
1169             rinv13           = gmx_mm_invsqrt_ps(rsq13);
1170             rinv21           = gmx_mm_invsqrt_ps(rsq21);
1171             rinv22           = gmx_mm_invsqrt_ps(rsq22);
1172             rinv23           = gmx_mm_invsqrt_ps(rsq23);
1173             rinv31           = gmx_mm_invsqrt_ps(rsq31);
1174             rinv32           = gmx_mm_invsqrt_ps(rsq32);
1175             rinv33           = gmx_mm_invsqrt_ps(rsq33);
1176
1177             rinvsq00         = gmx_mm_inv_ps(rsq00);
1178             rinvsq11         = _mm_mul_ps(rinv11,rinv11);
1179             rinvsq12         = _mm_mul_ps(rinv12,rinv12);
1180             rinvsq13         = _mm_mul_ps(rinv13,rinv13);
1181             rinvsq21         = _mm_mul_ps(rinv21,rinv21);
1182             rinvsq22         = _mm_mul_ps(rinv22,rinv22);
1183             rinvsq23         = _mm_mul_ps(rinv23,rinv23);
1184             rinvsq31         = _mm_mul_ps(rinv31,rinv31);
1185             rinvsq32         = _mm_mul_ps(rinv32,rinv32);
1186             rinvsq33         = _mm_mul_ps(rinv33,rinv33);
1187
1188             fjx0             = _mm_setzero_ps();
1189             fjy0             = _mm_setzero_ps();
1190             fjz0             = _mm_setzero_ps();
1191             fjx1             = _mm_setzero_ps();
1192             fjy1             = _mm_setzero_ps();
1193             fjz1             = _mm_setzero_ps();
1194             fjx2             = _mm_setzero_ps();
1195             fjy2             = _mm_setzero_ps();
1196             fjz2             = _mm_setzero_ps();
1197             fjx3             = _mm_setzero_ps();
1198             fjy3             = _mm_setzero_ps();
1199             fjz3             = _mm_setzero_ps();
1200
1201             /**************************
1202              * CALCULATE INTERACTIONS *
1203              **************************/
1204
1205             /* LENNARD-JONES DISPERSION/REPULSION */
1206
1207             rinvsix          = _mm_mul_ps(_mm_mul_ps(rinvsq00,rinvsq00),rinvsq00);
1208             fvdw             = _mm_mul_ps(_mm_msub_ps(c12_00,rinvsix,c6_00),_mm_mul_ps(rinvsix,rinvsq00));
1209
1210             fscal            = fvdw;
1211
1212              /* Update vectorial force */
1213             fix0             = _mm_macc_ps(dx00,fscal,fix0);
1214             fiy0             = _mm_macc_ps(dy00,fscal,fiy0);
1215             fiz0             = _mm_macc_ps(dz00,fscal,fiz0);
1216
1217             fjx0             = _mm_macc_ps(dx00,fscal,fjx0);
1218             fjy0             = _mm_macc_ps(dy00,fscal,fjy0);
1219             fjz0             = _mm_macc_ps(dz00,fscal,fjz0);
1220
1221             /**************************
1222              * CALCULATE INTERACTIONS *
1223              **************************/
1224
1225             /* REACTION-FIELD ELECTROSTATICS */
1226             felec            = _mm_mul_ps(qq11,_mm_msub_ps(rinv11,rinvsq11,krf2));
1227
1228             fscal            = felec;
1229
1230              /* Update vectorial force */
1231             fix1             = _mm_macc_ps(dx11,fscal,fix1);
1232             fiy1             = _mm_macc_ps(dy11,fscal,fiy1);
1233             fiz1             = _mm_macc_ps(dz11,fscal,fiz1);
1234
1235             fjx1             = _mm_macc_ps(dx11,fscal,fjx1);
1236             fjy1             = _mm_macc_ps(dy11,fscal,fjy1);
1237             fjz1             = _mm_macc_ps(dz11,fscal,fjz1);
1238
1239             /**************************
1240              * CALCULATE INTERACTIONS *
1241              **************************/
1242
1243             /* REACTION-FIELD ELECTROSTATICS */
1244             felec            = _mm_mul_ps(qq12,_mm_msub_ps(rinv12,rinvsq12,krf2));
1245
1246             fscal            = felec;
1247
1248              /* Update vectorial force */
1249             fix1             = _mm_macc_ps(dx12,fscal,fix1);
1250             fiy1             = _mm_macc_ps(dy12,fscal,fiy1);
1251             fiz1             = _mm_macc_ps(dz12,fscal,fiz1);
1252
1253             fjx2             = _mm_macc_ps(dx12,fscal,fjx2);
1254             fjy2             = _mm_macc_ps(dy12,fscal,fjy2);
1255             fjz2             = _mm_macc_ps(dz12,fscal,fjz2);
1256
1257             /**************************
1258              * CALCULATE INTERACTIONS *
1259              **************************/
1260
1261             /* REACTION-FIELD ELECTROSTATICS */
1262             felec            = _mm_mul_ps(qq13,_mm_msub_ps(rinv13,rinvsq13,krf2));
1263
1264             fscal            = felec;
1265
1266              /* Update vectorial force */
1267             fix1             = _mm_macc_ps(dx13,fscal,fix1);
1268             fiy1             = _mm_macc_ps(dy13,fscal,fiy1);
1269             fiz1             = _mm_macc_ps(dz13,fscal,fiz1);
1270
1271             fjx3             = _mm_macc_ps(dx13,fscal,fjx3);
1272             fjy3             = _mm_macc_ps(dy13,fscal,fjy3);
1273             fjz3             = _mm_macc_ps(dz13,fscal,fjz3);
1274
1275             /**************************
1276              * CALCULATE INTERACTIONS *
1277              **************************/
1278
1279             /* REACTION-FIELD ELECTROSTATICS */
1280             felec            = _mm_mul_ps(qq21,_mm_msub_ps(rinv21,rinvsq21,krf2));
1281
1282             fscal            = felec;
1283
1284              /* Update vectorial force */
1285             fix2             = _mm_macc_ps(dx21,fscal,fix2);
1286             fiy2             = _mm_macc_ps(dy21,fscal,fiy2);
1287             fiz2             = _mm_macc_ps(dz21,fscal,fiz2);
1288
1289             fjx1             = _mm_macc_ps(dx21,fscal,fjx1);
1290             fjy1             = _mm_macc_ps(dy21,fscal,fjy1);
1291             fjz1             = _mm_macc_ps(dz21,fscal,fjz1);
1292
1293             /**************************
1294              * CALCULATE INTERACTIONS *
1295              **************************/
1296
1297             /* REACTION-FIELD ELECTROSTATICS */
1298             felec            = _mm_mul_ps(qq22,_mm_msub_ps(rinv22,rinvsq22,krf2));
1299
1300             fscal            = felec;
1301
1302              /* Update vectorial force */
1303             fix2             = _mm_macc_ps(dx22,fscal,fix2);
1304             fiy2             = _mm_macc_ps(dy22,fscal,fiy2);
1305             fiz2             = _mm_macc_ps(dz22,fscal,fiz2);
1306
1307             fjx2             = _mm_macc_ps(dx22,fscal,fjx2);
1308             fjy2             = _mm_macc_ps(dy22,fscal,fjy2);
1309             fjz2             = _mm_macc_ps(dz22,fscal,fjz2);
1310
1311             /**************************
1312              * CALCULATE INTERACTIONS *
1313              **************************/
1314
1315             /* REACTION-FIELD ELECTROSTATICS */
1316             felec            = _mm_mul_ps(qq23,_mm_msub_ps(rinv23,rinvsq23,krf2));
1317
1318             fscal            = felec;
1319
1320              /* Update vectorial force */
1321             fix2             = _mm_macc_ps(dx23,fscal,fix2);
1322             fiy2             = _mm_macc_ps(dy23,fscal,fiy2);
1323             fiz2             = _mm_macc_ps(dz23,fscal,fiz2);
1324
1325             fjx3             = _mm_macc_ps(dx23,fscal,fjx3);
1326             fjy3             = _mm_macc_ps(dy23,fscal,fjy3);
1327             fjz3             = _mm_macc_ps(dz23,fscal,fjz3);
1328
1329             /**************************
1330              * CALCULATE INTERACTIONS *
1331              **************************/
1332
1333             /* REACTION-FIELD ELECTROSTATICS */
1334             felec            = _mm_mul_ps(qq31,_mm_msub_ps(rinv31,rinvsq31,krf2));
1335
1336             fscal            = felec;
1337
1338              /* Update vectorial force */
1339             fix3             = _mm_macc_ps(dx31,fscal,fix3);
1340             fiy3             = _mm_macc_ps(dy31,fscal,fiy3);
1341             fiz3             = _mm_macc_ps(dz31,fscal,fiz3);
1342
1343             fjx1             = _mm_macc_ps(dx31,fscal,fjx1);
1344             fjy1             = _mm_macc_ps(dy31,fscal,fjy1);
1345             fjz1             = _mm_macc_ps(dz31,fscal,fjz1);
1346
1347             /**************************
1348              * CALCULATE INTERACTIONS *
1349              **************************/
1350
1351             /* REACTION-FIELD ELECTROSTATICS */
1352             felec            = _mm_mul_ps(qq32,_mm_msub_ps(rinv32,rinvsq32,krf2));
1353
1354             fscal            = felec;
1355
1356              /* Update vectorial force */
1357             fix3             = _mm_macc_ps(dx32,fscal,fix3);
1358             fiy3             = _mm_macc_ps(dy32,fscal,fiy3);
1359             fiz3             = _mm_macc_ps(dz32,fscal,fiz3);
1360
1361             fjx2             = _mm_macc_ps(dx32,fscal,fjx2);
1362             fjy2             = _mm_macc_ps(dy32,fscal,fjy2);
1363             fjz2             = _mm_macc_ps(dz32,fscal,fjz2);
1364
1365             /**************************
1366              * CALCULATE INTERACTIONS *
1367              **************************/
1368
1369             /* REACTION-FIELD ELECTROSTATICS */
1370             felec            = _mm_mul_ps(qq33,_mm_msub_ps(rinv33,rinvsq33,krf2));
1371
1372             fscal            = felec;
1373
1374              /* Update vectorial force */
1375             fix3             = _mm_macc_ps(dx33,fscal,fix3);
1376             fiy3             = _mm_macc_ps(dy33,fscal,fiy3);
1377             fiz3             = _mm_macc_ps(dz33,fscal,fiz3);
1378
1379             fjx3             = _mm_macc_ps(dx33,fscal,fjx3);
1380             fjy3             = _mm_macc_ps(dy33,fscal,fjy3);
1381             fjz3             = _mm_macc_ps(dz33,fscal,fjz3);
1382
1383             fjptrA             = f+j_coord_offsetA;
1384             fjptrB             = f+j_coord_offsetB;
1385             fjptrC             = f+j_coord_offsetC;
1386             fjptrD             = f+j_coord_offsetD;
1387
1388             gmx_mm_decrement_4rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,
1389                                                    fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,
1390                                                    fjx2,fjy2,fjz2,fjx3,fjy3,fjz3);
1391
1392             /* Inner loop uses 303 flops */
1393         }
1394
1395         if(jidx<j_index_end)
1396         {
1397
1398             /* Get j neighbor index, and coordinate index */
1399             jnrlistA         = jjnr[jidx];
1400             jnrlistB         = jjnr[jidx+1];
1401             jnrlistC         = jjnr[jidx+2];
1402             jnrlistD         = jjnr[jidx+3];
1403             /* Sign of each element will be negative for non-real atoms.
1404              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
1405              * so use it as val = _mm_andnot_ps(mask,val) to clear dummy entries.
1406              */
1407             dummy_mask = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
1408             jnrA       = (jnrlistA>=0) ? jnrlistA : 0;
1409             jnrB       = (jnrlistB>=0) ? jnrlistB : 0;
1410             jnrC       = (jnrlistC>=0) ? jnrlistC : 0;
1411             jnrD       = (jnrlistD>=0) ? jnrlistD : 0;
1412             j_coord_offsetA  = DIM*jnrA;
1413             j_coord_offsetB  = DIM*jnrB;
1414             j_coord_offsetC  = DIM*jnrC;
1415             j_coord_offsetD  = DIM*jnrD;
1416
1417             /* load j atom coordinates */
1418             gmx_mm_load_4rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
1419                                               x+j_coord_offsetC,x+j_coord_offsetD,
1420                                               &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,
1421                                               &jy2,&jz2,&jx3,&jy3,&jz3);
1422
1423             /* Calculate displacement vector */
1424             dx00             = _mm_sub_ps(ix0,jx0);
1425             dy00             = _mm_sub_ps(iy0,jy0);
1426             dz00             = _mm_sub_ps(iz0,jz0);
1427             dx11             = _mm_sub_ps(ix1,jx1);
1428             dy11             = _mm_sub_ps(iy1,jy1);
1429             dz11             = _mm_sub_ps(iz1,jz1);
1430             dx12             = _mm_sub_ps(ix1,jx2);
1431             dy12             = _mm_sub_ps(iy1,jy2);
1432             dz12             = _mm_sub_ps(iz1,jz2);
1433             dx13             = _mm_sub_ps(ix1,jx3);
1434             dy13             = _mm_sub_ps(iy1,jy3);
1435             dz13             = _mm_sub_ps(iz1,jz3);
1436             dx21             = _mm_sub_ps(ix2,jx1);
1437             dy21             = _mm_sub_ps(iy2,jy1);
1438             dz21             = _mm_sub_ps(iz2,jz1);
1439             dx22             = _mm_sub_ps(ix2,jx2);
1440             dy22             = _mm_sub_ps(iy2,jy2);
1441             dz22             = _mm_sub_ps(iz2,jz2);
1442             dx23             = _mm_sub_ps(ix2,jx3);
1443             dy23             = _mm_sub_ps(iy2,jy3);
1444             dz23             = _mm_sub_ps(iz2,jz3);
1445             dx31             = _mm_sub_ps(ix3,jx1);
1446             dy31             = _mm_sub_ps(iy3,jy1);
1447             dz31             = _mm_sub_ps(iz3,jz1);
1448             dx32             = _mm_sub_ps(ix3,jx2);
1449             dy32             = _mm_sub_ps(iy3,jy2);
1450             dz32             = _mm_sub_ps(iz3,jz2);
1451             dx33             = _mm_sub_ps(ix3,jx3);
1452             dy33             = _mm_sub_ps(iy3,jy3);
1453             dz33             = _mm_sub_ps(iz3,jz3);
1454
1455             /* Calculate squared distance and things based on it */
1456             rsq00            = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
1457             rsq11            = gmx_mm_calc_rsq_ps(dx11,dy11,dz11);
1458             rsq12            = gmx_mm_calc_rsq_ps(dx12,dy12,dz12);
1459             rsq13            = gmx_mm_calc_rsq_ps(dx13,dy13,dz13);
1460             rsq21            = gmx_mm_calc_rsq_ps(dx21,dy21,dz21);
1461             rsq22            = gmx_mm_calc_rsq_ps(dx22,dy22,dz22);
1462             rsq23            = gmx_mm_calc_rsq_ps(dx23,dy23,dz23);
1463             rsq31            = gmx_mm_calc_rsq_ps(dx31,dy31,dz31);
1464             rsq32            = gmx_mm_calc_rsq_ps(dx32,dy32,dz32);
1465             rsq33            = gmx_mm_calc_rsq_ps(dx33,dy33,dz33);
1466
1467             rinv11           = gmx_mm_invsqrt_ps(rsq11);
1468             rinv12           = gmx_mm_invsqrt_ps(rsq12);
1469             rinv13           = gmx_mm_invsqrt_ps(rsq13);
1470             rinv21           = gmx_mm_invsqrt_ps(rsq21);
1471             rinv22           = gmx_mm_invsqrt_ps(rsq22);
1472             rinv23           = gmx_mm_invsqrt_ps(rsq23);
1473             rinv31           = gmx_mm_invsqrt_ps(rsq31);
1474             rinv32           = gmx_mm_invsqrt_ps(rsq32);
1475             rinv33           = gmx_mm_invsqrt_ps(rsq33);
1476
1477             rinvsq00         = gmx_mm_inv_ps(rsq00);
1478             rinvsq11         = _mm_mul_ps(rinv11,rinv11);
1479             rinvsq12         = _mm_mul_ps(rinv12,rinv12);
1480             rinvsq13         = _mm_mul_ps(rinv13,rinv13);
1481             rinvsq21         = _mm_mul_ps(rinv21,rinv21);
1482             rinvsq22         = _mm_mul_ps(rinv22,rinv22);
1483             rinvsq23         = _mm_mul_ps(rinv23,rinv23);
1484             rinvsq31         = _mm_mul_ps(rinv31,rinv31);
1485             rinvsq32         = _mm_mul_ps(rinv32,rinv32);
1486             rinvsq33         = _mm_mul_ps(rinv33,rinv33);
1487
1488             fjx0             = _mm_setzero_ps();
1489             fjy0             = _mm_setzero_ps();
1490             fjz0             = _mm_setzero_ps();
1491             fjx1             = _mm_setzero_ps();
1492             fjy1             = _mm_setzero_ps();
1493             fjz1             = _mm_setzero_ps();
1494             fjx2             = _mm_setzero_ps();
1495             fjy2             = _mm_setzero_ps();
1496             fjz2             = _mm_setzero_ps();
1497             fjx3             = _mm_setzero_ps();
1498             fjy3             = _mm_setzero_ps();
1499             fjz3             = _mm_setzero_ps();
1500
1501             /**************************
1502              * CALCULATE INTERACTIONS *
1503              **************************/
1504
1505             /* LENNARD-JONES DISPERSION/REPULSION */
1506
1507             rinvsix          = _mm_mul_ps(_mm_mul_ps(rinvsq00,rinvsq00),rinvsq00);
1508             fvdw             = _mm_mul_ps(_mm_msub_ps(c12_00,rinvsix,c6_00),_mm_mul_ps(rinvsix,rinvsq00));
1509
1510             fscal            = fvdw;
1511
1512             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1513
1514              /* Update vectorial force */
1515             fix0             = _mm_macc_ps(dx00,fscal,fix0);
1516             fiy0             = _mm_macc_ps(dy00,fscal,fiy0);
1517             fiz0             = _mm_macc_ps(dz00,fscal,fiz0);
1518
1519             fjx0             = _mm_macc_ps(dx00,fscal,fjx0);
1520             fjy0             = _mm_macc_ps(dy00,fscal,fjy0);
1521             fjz0             = _mm_macc_ps(dz00,fscal,fjz0);
1522
1523             /**************************
1524              * CALCULATE INTERACTIONS *
1525              **************************/
1526
1527             /* REACTION-FIELD ELECTROSTATICS */
1528             felec            = _mm_mul_ps(qq11,_mm_msub_ps(rinv11,rinvsq11,krf2));
1529
1530             fscal            = felec;
1531
1532             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1533
1534              /* Update vectorial force */
1535             fix1             = _mm_macc_ps(dx11,fscal,fix1);
1536             fiy1             = _mm_macc_ps(dy11,fscal,fiy1);
1537             fiz1             = _mm_macc_ps(dz11,fscal,fiz1);
1538
1539             fjx1             = _mm_macc_ps(dx11,fscal,fjx1);
1540             fjy1             = _mm_macc_ps(dy11,fscal,fjy1);
1541             fjz1             = _mm_macc_ps(dz11,fscal,fjz1);
1542
1543             /**************************
1544              * CALCULATE INTERACTIONS *
1545              **************************/
1546
1547             /* REACTION-FIELD ELECTROSTATICS */
1548             felec            = _mm_mul_ps(qq12,_mm_msub_ps(rinv12,rinvsq12,krf2));
1549
1550             fscal            = felec;
1551
1552             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1553
1554              /* Update vectorial force */
1555             fix1             = _mm_macc_ps(dx12,fscal,fix1);
1556             fiy1             = _mm_macc_ps(dy12,fscal,fiy1);
1557             fiz1             = _mm_macc_ps(dz12,fscal,fiz1);
1558
1559             fjx2             = _mm_macc_ps(dx12,fscal,fjx2);
1560             fjy2             = _mm_macc_ps(dy12,fscal,fjy2);
1561             fjz2             = _mm_macc_ps(dz12,fscal,fjz2);
1562
1563             /**************************
1564              * CALCULATE INTERACTIONS *
1565              **************************/
1566
1567             /* REACTION-FIELD ELECTROSTATICS */
1568             felec            = _mm_mul_ps(qq13,_mm_msub_ps(rinv13,rinvsq13,krf2));
1569
1570             fscal            = felec;
1571
1572             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1573
1574              /* Update vectorial force */
1575             fix1             = _mm_macc_ps(dx13,fscal,fix1);
1576             fiy1             = _mm_macc_ps(dy13,fscal,fiy1);
1577             fiz1             = _mm_macc_ps(dz13,fscal,fiz1);
1578
1579             fjx3             = _mm_macc_ps(dx13,fscal,fjx3);
1580             fjy3             = _mm_macc_ps(dy13,fscal,fjy3);
1581             fjz3             = _mm_macc_ps(dz13,fscal,fjz3);
1582
1583             /**************************
1584              * CALCULATE INTERACTIONS *
1585              **************************/
1586
1587             /* REACTION-FIELD ELECTROSTATICS */
1588             felec            = _mm_mul_ps(qq21,_mm_msub_ps(rinv21,rinvsq21,krf2));
1589
1590             fscal            = felec;
1591
1592             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1593
1594              /* Update vectorial force */
1595             fix2             = _mm_macc_ps(dx21,fscal,fix2);
1596             fiy2             = _mm_macc_ps(dy21,fscal,fiy2);
1597             fiz2             = _mm_macc_ps(dz21,fscal,fiz2);
1598
1599             fjx1             = _mm_macc_ps(dx21,fscal,fjx1);
1600             fjy1             = _mm_macc_ps(dy21,fscal,fjy1);
1601             fjz1             = _mm_macc_ps(dz21,fscal,fjz1);
1602
1603             /**************************
1604              * CALCULATE INTERACTIONS *
1605              **************************/
1606
1607             /* REACTION-FIELD ELECTROSTATICS */
1608             felec            = _mm_mul_ps(qq22,_mm_msub_ps(rinv22,rinvsq22,krf2));
1609
1610             fscal            = felec;
1611
1612             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1613
1614              /* Update vectorial force */
1615             fix2             = _mm_macc_ps(dx22,fscal,fix2);
1616             fiy2             = _mm_macc_ps(dy22,fscal,fiy2);
1617             fiz2             = _mm_macc_ps(dz22,fscal,fiz2);
1618
1619             fjx2             = _mm_macc_ps(dx22,fscal,fjx2);
1620             fjy2             = _mm_macc_ps(dy22,fscal,fjy2);
1621             fjz2             = _mm_macc_ps(dz22,fscal,fjz2);
1622
1623             /**************************
1624              * CALCULATE INTERACTIONS *
1625              **************************/
1626
1627             /* REACTION-FIELD ELECTROSTATICS */
1628             felec            = _mm_mul_ps(qq23,_mm_msub_ps(rinv23,rinvsq23,krf2));
1629
1630             fscal            = felec;
1631
1632             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1633
1634              /* Update vectorial force */
1635             fix2             = _mm_macc_ps(dx23,fscal,fix2);
1636             fiy2             = _mm_macc_ps(dy23,fscal,fiy2);
1637             fiz2             = _mm_macc_ps(dz23,fscal,fiz2);
1638
1639             fjx3             = _mm_macc_ps(dx23,fscal,fjx3);
1640             fjy3             = _mm_macc_ps(dy23,fscal,fjy3);
1641             fjz3             = _mm_macc_ps(dz23,fscal,fjz3);
1642
1643             /**************************
1644              * CALCULATE INTERACTIONS *
1645              **************************/
1646
1647             /* REACTION-FIELD ELECTROSTATICS */
1648             felec            = _mm_mul_ps(qq31,_mm_msub_ps(rinv31,rinvsq31,krf2));
1649
1650             fscal            = felec;
1651
1652             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1653
1654              /* Update vectorial force */
1655             fix3             = _mm_macc_ps(dx31,fscal,fix3);
1656             fiy3             = _mm_macc_ps(dy31,fscal,fiy3);
1657             fiz3             = _mm_macc_ps(dz31,fscal,fiz3);
1658
1659             fjx1             = _mm_macc_ps(dx31,fscal,fjx1);
1660             fjy1             = _mm_macc_ps(dy31,fscal,fjy1);
1661             fjz1             = _mm_macc_ps(dz31,fscal,fjz1);
1662
1663             /**************************
1664              * CALCULATE INTERACTIONS *
1665              **************************/
1666
1667             /* REACTION-FIELD ELECTROSTATICS */
1668             felec            = _mm_mul_ps(qq32,_mm_msub_ps(rinv32,rinvsq32,krf2));
1669
1670             fscal            = felec;
1671
1672             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1673
1674              /* Update vectorial force */
1675             fix3             = _mm_macc_ps(dx32,fscal,fix3);
1676             fiy3             = _mm_macc_ps(dy32,fscal,fiy3);
1677             fiz3             = _mm_macc_ps(dz32,fscal,fiz3);
1678
1679             fjx2             = _mm_macc_ps(dx32,fscal,fjx2);
1680             fjy2             = _mm_macc_ps(dy32,fscal,fjy2);
1681             fjz2             = _mm_macc_ps(dz32,fscal,fjz2);
1682
1683             /**************************
1684              * CALCULATE INTERACTIONS *
1685              **************************/
1686
1687             /* REACTION-FIELD ELECTROSTATICS */
1688             felec            = _mm_mul_ps(qq33,_mm_msub_ps(rinv33,rinvsq33,krf2));
1689
1690             fscal            = felec;
1691
1692             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1693
1694              /* Update vectorial force */
1695             fix3             = _mm_macc_ps(dx33,fscal,fix3);
1696             fiy3             = _mm_macc_ps(dy33,fscal,fiy3);
1697             fiz3             = _mm_macc_ps(dz33,fscal,fiz3);
1698
1699             fjx3             = _mm_macc_ps(dx33,fscal,fjx3);
1700             fjy3             = _mm_macc_ps(dy33,fscal,fjy3);
1701             fjz3             = _mm_macc_ps(dz33,fscal,fjz3);
1702
1703             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
1704             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
1705             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
1706             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
1707
1708             gmx_mm_decrement_4rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,
1709                                                    fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,
1710                                                    fjx2,fjy2,fjz2,fjx3,fjy3,fjz3);
1711
1712             /* Inner loop uses 303 flops */
1713         }
1714
1715         /* End of innermost loop */
1716
1717         gmx_mm_update_iforce_4atom_swizzle_ps(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
1718                                               f+i_coord_offset,fshift+i_shift_offset);
1719
1720         /* Increment number of inner iterations */
1721         inneriter                  += j_index_end - j_index_start;
1722
1723         /* Outer loop uses 24 flops */
1724     }
1725
1726     /* Increment number of outer iterations */
1727     outeriter        += nri;
1728
1729     /* Update outer/inner flops */
1730
1731     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4W4_F,outeriter*24 + inneriter*303);
1732 }