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