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

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