src/gromacs/gmxlib/nonbonded/nb_kernel_sse4_1_double/nb_kernel_ElecRFCut_VdwLJSw_GeomW4P1_sse4_1_double.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 sse4_1_double kernel generator.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "config.h"
  41
  42 #include <math.h>
  43
  44 #include "../nb_kernel.h"
  45 #include "gromacs/legacyheaders/types/simple.h"
  46 #include "gromacs/math/vec.h"
  47 #include "gromacs/legacyheaders/nrnb.h"
  48
  49 #include "gromacs/simd/math_x86_sse4_1_double.h"
  50 #include "kernelutil_x86_sse4_1_double.h"
  51
  52 /*
  53  * Gromacs nonbonded kernel:   nb_kernel_ElecRFCut_VdwLJSw_GeomW4P1_VF_sse4_1_double
  54  * Electrostatics interaction: ReactionField
  55  * VdW interaction:            LennardJones
  56  * Geometry:                   Water4-Particle
  57  * Calculate force/pot:        PotentialAndForce
  58  */
  59 void
  60 nb_kernel_ElecRFCut_VdwLJSw_GeomW4P1_VF_sse4_1_double
  61                     (t_nblist                    * gmx_restrict       nlist,
  62                      rvec                        * gmx_restrict          xx,
  63                      rvec                        * gmx_restrict          ff,
  64                      t_forcerec                  * gmx_restrict          fr,
  65                      t_mdatoms                   * gmx_restrict     mdatoms,
  66                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
  67                      t_nrnb                      * gmx_restrict        nrnb)
  68 {
  69     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
  70      * just 0 for non-waters.
  71      * Suffixes A,B refer to j loop unrolling done with SSE double precision, e.g. for the two different
  72      * jnr indices corresponding to data put in the four positions in the SIMD register.
  73      */
  74     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
  75     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
  76     int              jnrA,jnrB;
  77     int              j_coord_offsetA,j_coord_offsetB;
  78     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
  79     real             rcutoff_scalar;
  80     real             *shiftvec,*fshift,*x,*f;
  81     __m128d          tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
  82     int              vdwioffset0;
  83     __m128d          ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
  84     int              vdwioffset1;
  85     __m128d          ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
  86     int              vdwioffset2;
  87     __m128d          ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
  88     int              vdwioffset3;
  89     __m128d          ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
  90     int              vdwjidx0A,vdwjidx0B;
  91     __m128d          jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
  92     __m128d          dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
  93     __m128d          dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
  94     __m128d          dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
  95     __m128d          dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
  96     __m128d          velec,felec,velecsum,facel,crf,krf,krf2;
  97     real             *charge;
  98     int              nvdwtype;
  99     __m128d          rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
 100     int              *vdwtype;
 101     real             *vdwparam;
 102     __m128d          one_sixth   = _mm_set1_pd(1.0/6.0);
 103     __m128d          one_twelfth = _mm_set1_pd(1.0/12.0);
 104     __m128d          rswitch,swV3,swV4,swV5,swF2,swF3,swF4,d,d2,sw,dsw;
 105     real             rswitch_scalar,d_scalar;
 106     __m128d          dummy_mask,cutoff_mask;
 107     __m128d          signbit   = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
 108     __m128d          one     = _mm_set1_pd(1.0);
 109     __m128d          two     = _mm_set1_pd(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_pd(fr->epsfac);
 122     charge           = mdatoms->chargeA;
 123     krf              = _mm_set1_pd(fr->ic->k_rf);
 124     krf2             = _mm_set1_pd(fr->ic->k_rf*2.0);
 125     crf              = _mm_set1_pd(fr->ic->c_rf);
 126     nvdwtype         = fr->ntype;
 127     vdwparam         = fr->nbfp;
 128     vdwtype          = mdatoms->typeA;
 129
 130     /* Setup water-specific parameters */
 131     inr              = nlist->iinr[0];
 132     iq1              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+1]));
 133     iq2              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+2]));
 134     iq3              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+3]));
 135     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
 136
 137     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
 138     rcutoff_scalar   = fr->rcoulomb;
 139     rcutoff          = _mm_set1_pd(rcutoff_scalar);
 140     rcutoff2         = _mm_mul_pd(rcutoff,rcutoff);
 141
 142     rswitch_scalar   = fr->rvdw_switch;
 143     rswitch          = _mm_set1_pd(rswitch_scalar);
 144     /* Setup switch parameters */
 145     d_scalar         = rcutoff_scalar-rswitch_scalar;
 146     d                = _mm_set1_pd(d_scalar);
 147     swV3             = _mm_set1_pd(-10.0/(d_scalar*d_scalar*d_scalar));
 148     swV4             = _mm_set1_pd( 15.0/(d_scalar*d_scalar*d_scalar*d_scalar));
 149     swV5             = _mm_set1_pd( -6.0/(d_scalar*d_scalar*d_scalar*d_scalar*d_scalar));
 150     swF2             = _mm_set1_pd(-30.0/(d_scalar*d_scalar*d_scalar));
 151     swF3             = _mm_set1_pd( 60.0/(d_scalar*d_scalar*d_scalar*d_scalar));
 152     swF4             = _mm_set1_pd(-30.0/(d_scalar*d_scalar*d_scalar*d_scalar*d_scalar));
 153
 154     /* Avoid stupid compiler warnings */
 155     jnrA = jnrB = 0;
 156     j_coord_offsetA = 0;
 157     j_coord_offsetB = 0;
 158
 159     outeriter        = 0;
 160     inneriter        = 0;
 161
 162     /* Start outer loop over neighborlists */
 163     for(iidx=0; iidx<nri; iidx++)
 164     {
 165         /* Load shift vector for this list */
 166         i_shift_offset   = DIM*shiftidx[iidx];
 167
 168         /* Load limits for loop over neighbors */
 169         j_index_start    = jindex[iidx];
 170         j_index_end      = jindex[iidx+1];
 171
 172         /* Get outer coordinate index */
 173         inr              = iinr[iidx];
 174         i_coord_offset   = DIM*inr;
 175
 176         /* Load i particle coords and add shift vector */
 177         gmx_mm_load_shift_and_4rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,
 178                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
 179
 180         fix0             = _mm_setzero_pd();
 181         fiy0             = _mm_setzero_pd();
 182         fiz0             = _mm_setzero_pd();
 183         fix1             = _mm_setzero_pd();
 184         fiy1             = _mm_setzero_pd();
 185         fiz1             = _mm_setzero_pd();
 186         fix2             = _mm_setzero_pd();
 187         fiy2             = _mm_setzero_pd();
 188         fiz2             = _mm_setzero_pd();
 189         fix3             = _mm_setzero_pd();
 190         fiy3             = _mm_setzero_pd();
 191         fiz3             = _mm_setzero_pd();
 192
 193         /* Reset potential sums */
 194         velecsum         = _mm_setzero_pd();
 195         vvdwsum          = _mm_setzero_pd();
 196
 197         /* Start inner kernel loop */
 198         for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
 199         {
 200
 201             /* Get j neighbor index, and coordinate index */
 202             jnrA             = jjnr[jidx];
 203             jnrB             = jjnr[jidx+1];
 204             j_coord_offsetA  = DIM*jnrA;
 205             j_coord_offsetB  = DIM*jnrB;
 206
 207             /* load j atom coordinates */
 208             gmx_mm_load_1rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
 209                                               &jx0,&jy0,&jz0);
 210
 211             /* Calculate displacement vector */
 212             dx00             = _mm_sub_pd(ix0,jx0);
 213             dy00             = _mm_sub_pd(iy0,jy0);
 214             dz00             = _mm_sub_pd(iz0,jz0);
 215             dx10             = _mm_sub_pd(ix1,jx0);
 216             dy10             = _mm_sub_pd(iy1,jy0);
 217             dz10             = _mm_sub_pd(iz1,jz0);
 218             dx20             = _mm_sub_pd(ix2,jx0);
 219             dy20             = _mm_sub_pd(iy2,jy0);
 220             dz20             = _mm_sub_pd(iz2,jz0);
 221             dx30             = _mm_sub_pd(ix3,jx0);
 222             dy30             = _mm_sub_pd(iy3,jy0);
 223             dz30             = _mm_sub_pd(iz3,jz0);
 224
 225             /* Calculate squared distance and things based on it */
 226             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
 227             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
 228             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
 229             rsq30            = gmx_mm_calc_rsq_pd(dx30,dy30,dz30);
 230
 231             rinv00           = gmx_mm_invsqrt_pd(rsq00);
 232             rinv10           = gmx_mm_invsqrt_pd(rsq10);
 233             rinv20           = gmx_mm_invsqrt_pd(rsq20);
 234             rinv30           = gmx_mm_invsqrt_pd(rsq30);
 235
 236             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
 237             rinvsq10         = _mm_mul_pd(rinv10,rinv10);
 238             rinvsq20         = _mm_mul_pd(rinv20,rinv20);
 239             rinvsq30         = _mm_mul_pd(rinv30,rinv30);
 240
 241             /* Load parameters for j particles */
 242             jq0              = gmx_mm_load_2real_swizzle_pd(charge+jnrA+0,charge+jnrB+0);
 243             vdwjidx0A        = 2*vdwtype[jnrA+0];
 244             vdwjidx0B        = 2*vdwtype[jnrB+0];
 245
 246             fjx0             = _mm_setzero_pd();
 247             fjy0             = _mm_setzero_pd();
 248             fjz0             = _mm_setzero_pd();
 249
 250             /**************************
 251              * CALCULATE INTERACTIONS *
 252              **************************/
 253
 254             if (gmx_mm_any_lt(rsq00,rcutoff2))
 255             {
 256
 257             r00              = _mm_mul_pd(rsq00,rinv00);
 258
 259             /* Compute parameters for interactions between i and j atoms */
 260             gmx_mm_load_2pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,
 261                                          vdwparam+vdwioffset0+vdwjidx0B,&c6_00,&c12_00);
 262
 263             /* LENNARD-JONES DISPERSION/REPULSION */
 264
 265             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
 266             vvdw6            = _mm_mul_pd(c6_00,rinvsix);
 267             vvdw12           = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
 268             vvdw             = _mm_sub_pd( _mm_mul_pd(vvdw12,one_twelfth) , _mm_mul_pd(vvdw6,one_sixth) );
 269             fvdw             = _mm_mul_pd(_mm_sub_pd(vvdw12,vvdw6),rinvsq00);
 270
 271             d                = _mm_sub_pd(r00,rswitch);
 272             d                = _mm_max_pd(d,_mm_setzero_pd());
 273             d2               = _mm_mul_pd(d,d);
 274             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
 275
 276             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
 277
 278             /* Evaluate switch function */
 279             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 280             fvdw             = _mm_sub_pd( _mm_mul_pd(fvdw,sw) , _mm_mul_pd(rinv00,_mm_mul_pd(vvdw,dsw)) );
 281             vvdw             = _mm_mul_pd(vvdw,sw);
 282             cutoff_mask      = _mm_cmplt_pd(rsq00,rcutoff2);
 283
 284             /* Update potential sum for this i atom from the interaction with this j atom. */
 285             vvdw             = _mm_and_pd(vvdw,cutoff_mask);
 286             vvdwsum          = _mm_add_pd(vvdwsum,vvdw);
 287
 288             fscal            = fvdw;
 289
 290             fscal            = _mm_and_pd(fscal,cutoff_mask);
 291
 292             /* Calculate temporary vectorial force */
 293             tx               = _mm_mul_pd(fscal,dx00);
 294             ty               = _mm_mul_pd(fscal,dy00);
 295             tz               = _mm_mul_pd(fscal,dz00);
 296
 297             /* Update vectorial force */
 298             fix0             = _mm_add_pd(fix0,tx);
 299             fiy0             = _mm_add_pd(fiy0,ty);
 300             fiz0             = _mm_add_pd(fiz0,tz);
 301
 302             fjx0             = _mm_add_pd(fjx0,tx);
 303             fjy0             = _mm_add_pd(fjy0,ty);
 304             fjz0             = _mm_add_pd(fjz0,tz);
 305
 306             }
 307
 308             /**************************
 309              * CALCULATE INTERACTIONS *
 310              **************************/
 311
 312             if (gmx_mm_any_lt(rsq10,rcutoff2))
 313             {
 314
 315             /* Compute parameters for interactions between i and j atoms */
 316             qq10             = _mm_mul_pd(iq1,jq0);
 317
 318             /* REACTION-FIELD ELECTROSTATICS */
 319             velec            = _mm_mul_pd(qq10,_mm_sub_pd(_mm_add_pd(rinv10,_mm_mul_pd(krf,rsq10)),crf));
 320             felec            = _mm_mul_pd(qq10,_mm_sub_pd(_mm_mul_pd(rinv10,rinvsq10),krf2));
 321
 322             cutoff_mask      = _mm_cmplt_pd(rsq10,rcutoff2);
 323
 324             /* Update potential sum for this i atom from the interaction with this j atom. */
 325             velec            = _mm_and_pd(velec,cutoff_mask);
 326             velecsum         = _mm_add_pd(velecsum,velec);
 327
 328             fscal            = felec;
 329
 330             fscal            = _mm_and_pd(fscal,cutoff_mask);
 331
 332             /* Calculate temporary vectorial force */
 333             tx               = _mm_mul_pd(fscal,dx10);
 334             ty               = _mm_mul_pd(fscal,dy10);
 335             tz               = _mm_mul_pd(fscal,dz10);
 336
 337             /* Update vectorial force */
 338             fix1             = _mm_add_pd(fix1,tx);
 339             fiy1             = _mm_add_pd(fiy1,ty);
 340             fiz1             = _mm_add_pd(fiz1,tz);
 341
 342             fjx0             = _mm_add_pd(fjx0,tx);
 343             fjy0             = _mm_add_pd(fjy0,ty);
 344             fjz0             = _mm_add_pd(fjz0,tz);
 345
 346             }
 347
 348             /**************************
 349              * CALCULATE INTERACTIONS *
 350              **************************/
 351
 352             if (gmx_mm_any_lt(rsq20,rcutoff2))
 353             {
 354
 355             /* Compute parameters for interactions between i and j atoms */
 356             qq20             = _mm_mul_pd(iq2,jq0);
 357
 358             /* REACTION-FIELD ELECTROSTATICS */
 359             velec            = _mm_mul_pd(qq20,_mm_sub_pd(_mm_add_pd(rinv20,_mm_mul_pd(krf,rsq20)),crf));
 360             felec            = _mm_mul_pd(qq20,_mm_sub_pd(_mm_mul_pd(rinv20,rinvsq20),krf2));
 361
 362             cutoff_mask      = _mm_cmplt_pd(rsq20,rcutoff2);
 363
 364             /* Update potential sum for this i atom from the interaction with this j atom. */
 365             velec            = _mm_and_pd(velec,cutoff_mask);
 366             velecsum         = _mm_add_pd(velecsum,velec);
 367
 368             fscal            = felec;
 369
 370             fscal            = _mm_and_pd(fscal,cutoff_mask);
 371
 372             /* Calculate temporary vectorial force */
 373             tx               = _mm_mul_pd(fscal,dx20);
 374             ty               = _mm_mul_pd(fscal,dy20);
 375             tz               = _mm_mul_pd(fscal,dz20);
 376
 377             /* Update vectorial force */
 378             fix2             = _mm_add_pd(fix2,tx);
 379             fiy2             = _mm_add_pd(fiy2,ty);
 380             fiz2             = _mm_add_pd(fiz2,tz);
 381
 382             fjx0             = _mm_add_pd(fjx0,tx);
 383             fjy0             = _mm_add_pd(fjy0,ty);
 384             fjz0             = _mm_add_pd(fjz0,tz);
 385
 386             }
 387
 388             /**************************
 389              * CALCULATE INTERACTIONS *
 390              **************************/
 391
 392             if (gmx_mm_any_lt(rsq30,rcutoff2))
 393             {
 394
 395             /* Compute parameters for interactions between i and j atoms */
 396             qq30             = _mm_mul_pd(iq3,jq0);
 397
 398             /* REACTION-FIELD ELECTROSTATICS */
 399             velec            = _mm_mul_pd(qq30,_mm_sub_pd(_mm_add_pd(rinv30,_mm_mul_pd(krf,rsq30)),crf));
 400             felec            = _mm_mul_pd(qq30,_mm_sub_pd(_mm_mul_pd(rinv30,rinvsq30),krf2));
 401
 402             cutoff_mask      = _mm_cmplt_pd(rsq30,rcutoff2);
 403
 404             /* Update potential sum for this i atom from the interaction with this j atom. */
 405             velec            = _mm_and_pd(velec,cutoff_mask);
 406             velecsum         = _mm_add_pd(velecsum,velec);
 407
 408             fscal            = felec;
 409
 410             fscal            = _mm_and_pd(fscal,cutoff_mask);
 411
 412             /* Calculate temporary vectorial force */
 413             tx               = _mm_mul_pd(fscal,dx30);
 414             ty               = _mm_mul_pd(fscal,dy30);
 415             tz               = _mm_mul_pd(fscal,dz30);
 416
 417             /* Update vectorial force */
 418             fix3             = _mm_add_pd(fix3,tx);
 419             fiy3             = _mm_add_pd(fiy3,ty);
 420             fiz3             = _mm_add_pd(fiz3,tz);
 421
 422             fjx0             = _mm_add_pd(fjx0,tx);
 423             fjy0             = _mm_add_pd(fjy0,ty);
 424             fjz0             = _mm_add_pd(fjz0,tz);
 425
 426             }
 427
 428             gmx_mm_decrement_1rvec_2ptr_swizzle_pd(f+j_coord_offsetA,f+j_coord_offsetB,fjx0,fjy0,fjz0);
 429
 430             /* Inner loop uses 170 flops */
 431         }
 432
 433         if(jidx<j_index_end)
 434         {
 435
 436             jnrA             = jjnr[jidx];
 437             j_coord_offsetA  = DIM*jnrA;
 438
 439             /* load j atom coordinates */
 440             gmx_mm_load_1rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
 441                                               &jx0,&jy0,&jz0);
 442
 443             /* Calculate displacement vector */
 444             dx00             = _mm_sub_pd(ix0,jx0);
 445             dy00             = _mm_sub_pd(iy0,jy0);
 446             dz00             = _mm_sub_pd(iz0,jz0);
 447             dx10             = _mm_sub_pd(ix1,jx0);
 448             dy10             = _mm_sub_pd(iy1,jy0);
 449             dz10             = _mm_sub_pd(iz1,jz0);
 450             dx20             = _mm_sub_pd(ix2,jx0);
 451             dy20             = _mm_sub_pd(iy2,jy0);
 452             dz20             = _mm_sub_pd(iz2,jz0);
 453             dx30             = _mm_sub_pd(ix3,jx0);
 454             dy30             = _mm_sub_pd(iy3,jy0);
 455             dz30             = _mm_sub_pd(iz3,jz0);
 456
 457             /* Calculate squared distance and things based on it */
 458             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
 459             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
 460             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
 461             rsq30            = gmx_mm_calc_rsq_pd(dx30,dy30,dz30);
 462
 463             rinv00           = gmx_mm_invsqrt_pd(rsq00);
 464             rinv10           = gmx_mm_invsqrt_pd(rsq10);
 465             rinv20           = gmx_mm_invsqrt_pd(rsq20);
 466             rinv30           = gmx_mm_invsqrt_pd(rsq30);
 467
 468             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
 469             rinvsq10         = _mm_mul_pd(rinv10,rinv10);
 470             rinvsq20         = _mm_mul_pd(rinv20,rinv20);
 471             rinvsq30         = _mm_mul_pd(rinv30,rinv30);
 472
 473             /* Load parameters for j particles */
 474             jq0              = _mm_load_sd(charge+jnrA+0);
 475             vdwjidx0A        = 2*vdwtype[jnrA+0];
 476
 477             fjx0             = _mm_setzero_pd();
 478             fjy0             = _mm_setzero_pd();
 479             fjz0             = _mm_setzero_pd();
 480
 481             /**************************
 482              * CALCULATE INTERACTIONS *
 483              **************************/
 484
 485             if (gmx_mm_any_lt(rsq00,rcutoff2))
 486             {
 487
 488             r00              = _mm_mul_pd(rsq00,rinv00);
 489
 490             /* Compute parameters for interactions between i and j atoms */
 491             gmx_mm_load_1pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,&c6_00,&c12_00);
 492
 493             /* LENNARD-JONES DISPERSION/REPULSION */
 494
 495             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
 496             vvdw6            = _mm_mul_pd(c6_00,rinvsix);
 497             vvdw12           = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
 498             vvdw             = _mm_sub_pd( _mm_mul_pd(vvdw12,one_twelfth) , _mm_mul_pd(vvdw6,one_sixth) );
 499             fvdw             = _mm_mul_pd(_mm_sub_pd(vvdw12,vvdw6),rinvsq00);
 500
 501             d                = _mm_sub_pd(r00,rswitch);
 502             d                = _mm_max_pd(d,_mm_setzero_pd());
 503             d2               = _mm_mul_pd(d,d);
 504             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
 505
 506             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
 507
 508             /* Evaluate switch function */
 509             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 510             fvdw             = _mm_sub_pd( _mm_mul_pd(fvdw,sw) , _mm_mul_pd(rinv00,_mm_mul_pd(vvdw,dsw)) );
 511             vvdw             = _mm_mul_pd(vvdw,sw);
 512             cutoff_mask      = _mm_cmplt_pd(rsq00,rcutoff2);
 513
 514             /* Update potential sum for this i atom from the interaction with this j atom. */
 515             vvdw             = _mm_and_pd(vvdw,cutoff_mask);
 516             vvdw             = _mm_unpacklo_pd(vvdw,_mm_setzero_pd());
 517             vvdwsum          = _mm_add_pd(vvdwsum,vvdw);
 518
 519             fscal            = fvdw;
 520
 521             fscal            = _mm_and_pd(fscal,cutoff_mask);
 522
 523             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
 524
 525             /* Calculate temporary vectorial force */
 526             tx               = _mm_mul_pd(fscal,dx00);
 527             ty               = _mm_mul_pd(fscal,dy00);
 528             tz               = _mm_mul_pd(fscal,dz00);
 529
 530             /* Update vectorial force */
 531             fix0             = _mm_add_pd(fix0,tx);
 532             fiy0             = _mm_add_pd(fiy0,ty);
 533             fiz0             = _mm_add_pd(fiz0,tz);
 534
 535             fjx0             = _mm_add_pd(fjx0,tx);
 536             fjy0             = _mm_add_pd(fjy0,ty);
 537             fjz0             = _mm_add_pd(fjz0,tz);
 538
 539             }
 540
 541             /**************************
 542              * CALCULATE INTERACTIONS *
 543              **************************/
 544
 545             if (gmx_mm_any_lt(rsq10,rcutoff2))
 546             {
 547
 548             /* Compute parameters for interactions between i and j atoms */
 549             qq10             = _mm_mul_pd(iq1,jq0);
 550
 551             /* REACTION-FIELD ELECTROSTATICS */
 552             velec            = _mm_mul_pd(qq10,_mm_sub_pd(_mm_add_pd(rinv10,_mm_mul_pd(krf,rsq10)),crf));
 553             felec            = _mm_mul_pd(qq10,_mm_sub_pd(_mm_mul_pd(rinv10,rinvsq10),krf2));
 554
 555             cutoff_mask      = _mm_cmplt_pd(rsq10,rcutoff2);
 556
 557             /* Update potential sum for this i atom from the interaction with this j atom. */
 558             velec            = _mm_and_pd(velec,cutoff_mask);
 559             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
 560             velecsum         = _mm_add_pd(velecsum,velec);
 561
 562             fscal            = felec;
 563
 564             fscal            = _mm_and_pd(fscal,cutoff_mask);
 565
 566             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
 567
 568             /* Calculate temporary vectorial force */
 569             tx               = _mm_mul_pd(fscal,dx10);
 570             ty               = _mm_mul_pd(fscal,dy10);
 571             tz               = _mm_mul_pd(fscal,dz10);
 572
 573             /* Update vectorial force */
 574             fix1             = _mm_add_pd(fix1,tx);
 575             fiy1             = _mm_add_pd(fiy1,ty);
 576             fiz1             = _mm_add_pd(fiz1,tz);
 577
 578             fjx0             = _mm_add_pd(fjx0,tx);
 579             fjy0             = _mm_add_pd(fjy0,ty);
 580             fjz0             = _mm_add_pd(fjz0,tz);
 581
 582             }
 583
 584             /**************************
 585              * CALCULATE INTERACTIONS *
 586              **************************/
 587
 588             if (gmx_mm_any_lt(rsq20,rcutoff2))
 589             {
 590
 591             /* Compute parameters for interactions between i and j atoms */
 592             qq20             = _mm_mul_pd(iq2,jq0);
 593
 594             /* REACTION-FIELD ELECTROSTATICS */
 595             velec            = _mm_mul_pd(qq20,_mm_sub_pd(_mm_add_pd(rinv20,_mm_mul_pd(krf,rsq20)),crf));
 596             felec            = _mm_mul_pd(qq20,_mm_sub_pd(_mm_mul_pd(rinv20,rinvsq20),krf2));
 597
 598             cutoff_mask      = _mm_cmplt_pd(rsq20,rcutoff2);
 599
 600             /* Update potential sum for this i atom from the interaction with this j atom. */
 601             velec            = _mm_and_pd(velec,cutoff_mask);
 602             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
 603             velecsum         = _mm_add_pd(velecsum,velec);
 604
 605             fscal            = felec;
 606
 607             fscal            = _mm_and_pd(fscal,cutoff_mask);
 608
 609             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
 610
 611             /* Calculate temporary vectorial force */
 612             tx               = _mm_mul_pd(fscal,dx20);
 613             ty               = _mm_mul_pd(fscal,dy20);
 614             tz               = _mm_mul_pd(fscal,dz20);
 615
 616             /* Update vectorial force */
 617             fix2             = _mm_add_pd(fix2,tx);
 618             fiy2             = _mm_add_pd(fiy2,ty);
 619             fiz2             = _mm_add_pd(fiz2,tz);
 620
 621             fjx0             = _mm_add_pd(fjx0,tx);
 622             fjy0             = _mm_add_pd(fjy0,ty);
 623             fjz0             = _mm_add_pd(fjz0,tz);
 624
 625             }
 626
 627             /**************************
 628              * CALCULATE INTERACTIONS *
 629              **************************/
 630
 631             if (gmx_mm_any_lt(rsq30,rcutoff2))
 632             {
 633
 634             /* Compute parameters for interactions between i and j atoms */
 635             qq30             = _mm_mul_pd(iq3,jq0);
 636
 637             /* REACTION-FIELD ELECTROSTATICS */
 638             velec            = _mm_mul_pd(qq30,_mm_sub_pd(_mm_add_pd(rinv30,_mm_mul_pd(krf,rsq30)),crf));
 639             felec            = _mm_mul_pd(qq30,_mm_sub_pd(_mm_mul_pd(rinv30,rinvsq30),krf2));
 640
 641             cutoff_mask      = _mm_cmplt_pd(rsq30,rcutoff2);
 642
 643             /* Update potential sum for this i atom from the interaction with this j atom. */
 644             velec            = _mm_and_pd(velec,cutoff_mask);
 645             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
 646             velecsum         = _mm_add_pd(velecsum,velec);
 647
 648             fscal            = felec;
 649
 650             fscal            = _mm_and_pd(fscal,cutoff_mask);
 651
 652             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
 653
 654             /* Calculate temporary vectorial force */
 655             tx               = _mm_mul_pd(fscal,dx30);
 656             ty               = _mm_mul_pd(fscal,dy30);
 657             tz               = _mm_mul_pd(fscal,dz30);
 658
 659             /* Update vectorial force */
 660             fix3             = _mm_add_pd(fix3,tx);
 661             fiy3             = _mm_add_pd(fiy3,ty);
 662             fiz3             = _mm_add_pd(fiz3,tz);
 663
 664             fjx0             = _mm_add_pd(fjx0,tx);
 665             fjy0             = _mm_add_pd(fjy0,ty);
 666             fjz0             = _mm_add_pd(fjz0,tz);
 667
 668             }
 669
 670             gmx_mm_decrement_1rvec_1ptr_swizzle_pd(f+j_coord_offsetA,fjx0,fjy0,fjz0);
 671
 672             /* Inner loop uses 170 flops */
 673         }
 674
 675         /* End of innermost loop */
 676
 677         gmx_mm_update_iforce_4atom_swizzle_pd(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
 678                                               f+i_coord_offset,fshift+i_shift_offset);
 679
 680         ggid                        = gid[iidx];
 681         /* Update potential energies */
 682         gmx_mm_update_1pot_pd(velecsum,kernel_data->energygrp_elec+ggid);
 683         gmx_mm_update_1pot_pd(vvdwsum,kernel_data->energygrp_vdw+ggid);
 684
 685         /* Increment number of inner iterations */
 686         inneriter                  += j_index_end - j_index_start;
 687
 688         /* Outer loop uses 26 flops */
 689     }
 690
 691     /* Increment number of outer iterations */
 692     outeriter        += nri;
 693
 694     /* Update outer/inner flops */
 695
 696     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4_VF,outeriter*26 + inneriter*170);
 697 }
 698 /*
 699  * Gromacs nonbonded kernel:   nb_kernel_ElecRFCut_VdwLJSw_GeomW4P1_F_sse4_1_double
 700  * Electrostatics interaction: ReactionField
 701  * VdW interaction:            LennardJones
 702  * Geometry:                   Water4-Particle
 703  * Calculate force/pot:        Force
 704  */
 705 void
 706 nb_kernel_ElecRFCut_VdwLJSw_GeomW4P1_F_sse4_1_double
 707                     (t_nblist                    * gmx_restrict       nlist,
 708                      rvec                        * gmx_restrict          xx,
 709                      rvec                        * gmx_restrict          ff,
 710                      t_forcerec                  * gmx_restrict          fr,
 711                      t_mdatoms                   * gmx_restrict     mdatoms,
 712                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
 713                      t_nrnb                      * gmx_restrict        nrnb)
 714 {
 715     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
 716      * just 0 for non-waters.
 717      * Suffixes A,B refer to j loop unrolling done with SSE double precision, e.g. for the two different
 718      * jnr indices corresponding to data put in the four positions in the SIMD register.
 719      */
 720     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
 721     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
 722     int              jnrA,jnrB;
 723     int              j_coord_offsetA,j_coord_offsetB;
 724     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
 725     real             rcutoff_scalar;
 726     real             *shiftvec,*fshift,*x,*f;
 727     __m128d          tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
 728     int              vdwioffset0;
 729     __m128d          ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
 730     int              vdwioffset1;
 731     __m128d          ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
 732     int              vdwioffset2;
 733     __m128d          ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
 734     int              vdwioffset3;
 735     __m128d          ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
 736     int              vdwjidx0A,vdwjidx0B;
 737     __m128d          jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
 738     __m128d          dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
 739     __m128d          dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
 740     __m128d          dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
 741     __m128d          dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
 742     __m128d          velec,felec,velecsum,facel,crf,krf,krf2;
 743     real             *charge;
 744     int              nvdwtype;
 745     __m128d          rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
 746     int              *vdwtype;
 747     real             *vdwparam;
 748     __m128d          one_sixth   = _mm_set1_pd(1.0/6.0);
 749     __m128d          one_twelfth = _mm_set1_pd(1.0/12.0);
 750     __m128d          rswitch,swV3,swV4,swV5,swF2,swF3,swF4,d,d2,sw,dsw;
 751     real             rswitch_scalar,d_scalar;
 752     __m128d          dummy_mask,cutoff_mask;
 753     __m128d          signbit   = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
 754     __m128d          one     = _mm_set1_pd(1.0);
 755     __m128d          two     = _mm_set1_pd(2.0);
 756     x                = xx[0];
 757     f                = ff[0];
 758
 759     nri              = nlist->nri;
 760     iinr             = nlist->iinr;
 761     jindex           = nlist->jindex;
 762     jjnr             = nlist->jjnr;
 763     shiftidx         = nlist->shift;
 764     gid              = nlist->gid;
 765     shiftvec         = fr->shift_vec[0];
 766     fshift           = fr->fshift[0];
 767     facel            = _mm_set1_pd(fr->epsfac);
 768     charge           = mdatoms->chargeA;
 769     krf              = _mm_set1_pd(fr->ic->k_rf);
 770     krf2             = _mm_set1_pd(fr->ic->k_rf*2.0);
 771     crf              = _mm_set1_pd(fr->ic->c_rf);
 772     nvdwtype         = fr->ntype;
 773     vdwparam         = fr->nbfp;
 774     vdwtype          = mdatoms->typeA;
 775
 776     /* Setup water-specific parameters */
 777     inr              = nlist->iinr[0];
 778     iq1              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+1]));
 779     iq2              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+2]));
 780     iq3              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+3]));
 781     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
 782
 783     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
 784     rcutoff_scalar   = fr->rcoulomb;
 785     rcutoff          = _mm_set1_pd(rcutoff_scalar);
 786     rcutoff2         = _mm_mul_pd(rcutoff,rcutoff);
 787
 788     rswitch_scalar   = fr->rvdw_switch;
 789     rswitch          = _mm_set1_pd(rswitch_scalar);
 790     /* Setup switch parameters */
 791     d_scalar         = rcutoff_scalar-rswitch_scalar;
 792     d                = _mm_set1_pd(d_scalar);
 793     swV3             = _mm_set1_pd(-10.0/(d_scalar*d_scalar*d_scalar));
 794     swV4             = _mm_set1_pd( 15.0/(d_scalar*d_scalar*d_scalar*d_scalar));
 795     swV5             = _mm_set1_pd( -6.0/(d_scalar*d_scalar*d_scalar*d_scalar*d_scalar));
 796     swF2             = _mm_set1_pd(-30.0/(d_scalar*d_scalar*d_scalar));
 797     swF3             = _mm_set1_pd( 60.0/(d_scalar*d_scalar*d_scalar*d_scalar));
 798     swF4             = _mm_set1_pd(-30.0/(d_scalar*d_scalar*d_scalar*d_scalar*d_scalar));
 799
 800     /* Avoid stupid compiler warnings */
 801     jnrA = jnrB = 0;
 802     j_coord_offsetA = 0;
 803     j_coord_offsetB = 0;
 804
 805     outeriter        = 0;
 806     inneriter        = 0;
 807
 808     /* Start outer loop over neighborlists */
 809     for(iidx=0; iidx<nri; iidx++)
 810     {
 811         /* Load shift vector for this list */
 812         i_shift_offset   = DIM*shiftidx[iidx];
 813
 814         /* Load limits for loop over neighbors */
 815         j_index_start    = jindex[iidx];
 816         j_index_end      = jindex[iidx+1];
 817
 818         /* Get outer coordinate index */
 819         inr              = iinr[iidx];
 820         i_coord_offset   = DIM*inr;
 821
 822         /* Load i particle coords and add shift vector */
 823         gmx_mm_load_shift_and_4rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,
 824                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
 825
 826         fix0             = _mm_setzero_pd();
 827         fiy0             = _mm_setzero_pd();
 828         fiz0             = _mm_setzero_pd();
 829         fix1             = _mm_setzero_pd();
 830         fiy1             = _mm_setzero_pd();
 831         fiz1             = _mm_setzero_pd();
 832         fix2             = _mm_setzero_pd();
 833         fiy2             = _mm_setzero_pd();
 834         fiz2             = _mm_setzero_pd();
 835         fix3             = _mm_setzero_pd();
 836         fiy3             = _mm_setzero_pd();
 837         fiz3             = _mm_setzero_pd();
 838
 839         /* Start inner kernel loop */
 840         for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
 841         {
 842
 843             /* Get j neighbor index, and coordinate index */
 844             jnrA             = jjnr[jidx];
 845             jnrB             = jjnr[jidx+1];
 846             j_coord_offsetA  = DIM*jnrA;
 847             j_coord_offsetB  = DIM*jnrB;
 848
 849             /* load j atom coordinates */
 850             gmx_mm_load_1rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
 851                                               &jx0,&jy0,&jz0);
 852
 853             /* Calculate displacement vector */
 854             dx00             = _mm_sub_pd(ix0,jx0);
 855             dy00             = _mm_sub_pd(iy0,jy0);
 856             dz00             = _mm_sub_pd(iz0,jz0);
 857             dx10             = _mm_sub_pd(ix1,jx0);
 858             dy10             = _mm_sub_pd(iy1,jy0);
 859             dz10             = _mm_sub_pd(iz1,jz0);
 860             dx20             = _mm_sub_pd(ix2,jx0);
 861             dy20             = _mm_sub_pd(iy2,jy0);
 862             dz20             = _mm_sub_pd(iz2,jz0);
 863             dx30             = _mm_sub_pd(ix3,jx0);
 864             dy30             = _mm_sub_pd(iy3,jy0);
 865             dz30             = _mm_sub_pd(iz3,jz0);
 866
 867             /* Calculate squared distance and things based on it */
 868             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
 869             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
 870             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
 871             rsq30            = gmx_mm_calc_rsq_pd(dx30,dy30,dz30);
 872
 873             rinv00           = gmx_mm_invsqrt_pd(rsq00);
 874             rinv10           = gmx_mm_invsqrt_pd(rsq10);
 875             rinv20           = gmx_mm_invsqrt_pd(rsq20);
 876             rinv30           = gmx_mm_invsqrt_pd(rsq30);
 877
 878             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
 879             rinvsq10         = _mm_mul_pd(rinv10,rinv10);
 880             rinvsq20         = _mm_mul_pd(rinv20,rinv20);
 881             rinvsq30         = _mm_mul_pd(rinv30,rinv30);
 882
 883             /* Load parameters for j particles */
 884             jq0              = gmx_mm_load_2real_swizzle_pd(charge+jnrA+0,charge+jnrB+0);
 885             vdwjidx0A        = 2*vdwtype[jnrA+0];
 886             vdwjidx0B        = 2*vdwtype[jnrB+0];
 887
 888             fjx0             = _mm_setzero_pd();
 889             fjy0             = _mm_setzero_pd();
 890             fjz0             = _mm_setzero_pd();
 891
 892             /**************************
 893              * CALCULATE INTERACTIONS *
 894              **************************/
 895
 896             if (gmx_mm_any_lt(rsq00,rcutoff2))
 897             {
 898
 899             r00              = _mm_mul_pd(rsq00,rinv00);
 900
 901             /* Compute parameters for interactions between i and j atoms */
 902             gmx_mm_load_2pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,
 903                                          vdwparam+vdwioffset0+vdwjidx0B,&c6_00,&c12_00);
 904
 905             /* LENNARD-JONES DISPERSION/REPULSION */
 906
 907             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
 908             vvdw6            = _mm_mul_pd(c6_00,rinvsix);
 909             vvdw12           = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
 910             vvdw             = _mm_sub_pd( _mm_mul_pd(vvdw12,one_twelfth) , _mm_mul_pd(vvdw6,one_sixth) );
 911             fvdw             = _mm_mul_pd(_mm_sub_pd(vvdw12,vvdw6),rinvsq00);
 912
 913             d                = _mm_sub_pd(r00,rswitch);
 914             d                = _mm_max_pd(d,_mm_setzero_pd());
 915             d2               = _mm_mul_pd(d,d);
 916             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
 917
 918             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
 919
 920             /* Evaluate switch function */
 921             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
 922             fvdw             = _mm_sub_pd( _mm_mul_pd(fvdw,sw) , _mm_mul_pd(rinv00,_mm_mul_pd(vvdw,dsw)) );
 923             cutoff_mask      = _mm_cmplt_pd(rsq00,rcutoff2);
 924
 925             fscal            = fvdw;
 926
 927             fscal            = _mm_and_pd(fscal,cutoff_mask);
 928
 929             /* Calculate temporary vectorial force */
 930             tx               = _mm_mul_pd(fscal,dx00);
 931             ty               = _mm_mul_pd(fscal,dy00);
 932             tz               = _mm_mul_pd(fscal,dz00);
 933
 934             /* Update vectorial force */
 935             fix0             = _mm_add_pd(fix0,tx);
 936             fiy0             = _mm_add_pd(fiy0,ty);
 937             fiz0             = _mm_add_pd(fiz0,tz);
 938
 939             fjx0             = _mm_add_pd(fjx0,tx);
 940             fjy0             = _mm_add_pd(fjy0,ty);
 941             fjz0             = _mm_add_pd(fjz0,tz);
 942
 943             }
 944
 945             /**************************
 946              * CALCULATE INTERACTIONS *
 947              **************************/
 948
 949             if (gmx_mm_any_lt(rsq10,rcutoff2))
 950             {
 951
 952             /* Compute parameters for interactions between i and j atoms */
 953             qq10             = _mm_mul_pd(iq1,jq0);
 954
 955             /* REACTION-FIELD ELECTROSTATICS */
 956             felec            = _mm_mul_pd(qq10,_mm_sub_pd(_mm_mul_pd(rinv10,rinvsq10),krf2));
 957
 958             cutoff_mask      = _mm_cmplt_pd(rsq10,rcutoff2);
 959
 960             fscal            = felec;
 961
 962             fscal            = _mm_and_pd(fscal,cutoff_mask);
 963
 964             /* Calculate temporary vectorial force */
 965             tx               = _mm_mul_pd(fscal,dx10);
 966             ty               = _mm_mul_pd(fscal,dy10);
 967             tz               = _mm_mul_pd(fscal,dz10);
 968
 969             /* Update vectorial force */
 970             fix1             = _mm_add_pd(fix1,tx);
 971             fiy1             = _mm_add_pd(fiy1,ty);
 972             fiz1             = _mm_add_pd(fiz1,tz);
 973
 974             fjx0             = _mm_add_pd(fjx0,tx);
 975             fjy0             = _mm_add_pd(fjy0,ty);
 976             fjz0             = _mm_add_pd(fjz0,tz);
 977
 978             }
 979
 980             /**************************
 981              * CALCULATE INTERACTIONS *
 982              **************************/
 983
 984             if (gmx_mm_any_lt(rsq20,rcutoff2))
 985             {
 986
 987             /* Compute parameters for interactions between i and j atoms */
 988             qq20             = _mm_mul_pd(iq2,jq0);
 989
 990             /* REACTION-FIELD ELECTROSTATICS */
 991             felec            = _mm_mul_pd(qq20,_mm_sub_pd(_mm_mul_pd(rinv20,rinvsq20),krf2));
 992
 993             cutoff_mask      = _mm_cmplt_pd(rsq20,rcutoff2);
 994
 995             fscal            = felec;
 996
 997             fscal            = _mm_and_pd(fscal,cutoff_mask);
 998
 999             /* Calculate temporary vectorial force */
1000             tx               = _mm_mul_pd(fscal,dx20);
1001             ty               = _mm_mul_pd(fscal,dy20);
1002             tz               = _mm_mul_pd(fscal,dz20);
1003
1004             /* Update vectorial force */
1005             fix2             = _mm_add_pd(fix2,tx);
1006             fiy2             = _mm_add_pd(fiy2,ty);
1007             fiz2             = _mm_add_pd(fiz2,tz);
1008
1009             fjx0             = _mm_add_pd(fjx0,tx);
1010             fjy0             = _mm_add_pd(fjy0,ty);
1011             fjz0             = _mm_add_pd(fjz0,tz);
1012
1013             }
1014
1015             /**************************
1016              * CALCULATE INTERACTIONS *
1017              **************************/
1018
1019             if (gmx_mm_any_lt(rsq30,rcutoff2))
1020             {
1021
1022             /* Compute parameters for interactions between i and j atoms */
1023             qq30             = _mm_mul_pd(iq3,jq0);
1024
1025             /* REACTION-FIELD ELECTROSTATICS */
1026             felec            = _mm_mul_pd(qq30,_mm_sub_pd(_mm_mul_pd(rinv30,rinvsq30),krf2));
1027
1028             cutoff_mask      = _mm_cmplt_pd(rsq30,rcutoff2);
1029
1030             fscal            = felec;
1031
1032             fscal            = _mm_and_pd(fscal,cutoff_mask);
1033
1034             /* Calculate temporary vectorial force */
1035             tx               = _mm_mul_pd(fscal,dx30);
1036             ty               = _mm_mul_pd(fscal,dy30);
1037             tz               = _mm_mul_pd(fscal,dz30);
1038
1039             /* Update vectorial force */
1040             fix3             = _mm_add_pd(fix3,tx);
1041             fiy3             = _mm_add_pd(fiy3,ty);
1042             fiz3             = _mm_add_pd(fiz3,tz);
1043
1044             fjx0             = _mm_add_pd(fjx0,tx);
1045             fjy0             = _mm_add_pd(fjy0,ty);
1046             fjz0             = _mm_add_pd(fjz0,tz);
1047
1048             }
1049
1050             gmx_mm_decrement_1rvec_2ptr_swizzle_pd(f+j_coord_offsetA,f+j_coord_offsetB,fjx0,fjy0,fjz0);
1051
1052             /* Inner loop uses 149 flops */
1053         }
1054
1055         if(jidx<j_index_end)
1056         {
1057
1058             jnrA             = jjnr[jidx];
1059             j_coord_offsetA  = DIM*jnrA;
1060
1061             /* load j atom coordinates */
1062             gmx_mm_load_1rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
1063                                               &jx0,&jy0,&jz0);
1064
1065             /* Calculate displacement vector */
1066             dx00             = _mm_sub_pd(ix0,jx0);
1067             dy00             = _mm_sub_pd(iy0,jy0);
1068             dz00             = _mm_sub_pd(iz0,jz0);
1069             dx10             = _mm_sub_pd(ix1,jx0);
1070             dy10             = _mm_sub_pd(iy1,jy0);
1071             dz10             = _mm_sub_pd(iz1,jz0);
1072             dx20             = _mm_sub_pd(ix2,jx0);
1073             dy20             = _mm_sub_pd(iy2,jy0);
1074             dz20             = _mm_sub_pd(iz2,jz0);
1075             dx30             = _mm_sub_pd(ix3,jx0);
1076             dy30             = _mm_sub_pd(iy3,jy0);
1077             dz30             = _mm_sub_pd(iz3,jz0);
1078
1079             /* Calculate squared distance and things based on it */
1080             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
1081             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
1082             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
1083             rsq30            = gmx_mm_calc_rsq_pd(dx30,dy30,dz30);
1084
1085             rinv00           = gmx_mm_invsqrt_pd(rsq00);
1086             rinv10           = gmx_mm_invsqrt_pd(rsq10);
1087             rinv20           = gmx_mm_invsqrt_pd(rsq20);
1088             rinv30           = gmx_mm_invsqrt_pd(rsq30);
1089
1090             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
1091             rinvsq10         = _mm_mul_pd(rinv10,rinv10);
1092             rinvsq20         = _mm_mul_pd(rinv20,rinv20);
1093             rinvsq30         = _mm_mul_pd(rinv30,rinv30);
1094
1095             /* Load parameters for j particles */
1096             jq0              = _mm_load_sd(charge+jnrA+0);
1097             vdwjidx0A        = 2*vdwtype[jnrA+0];
1098
1099             fjx0             = _mm_setzero_pd();
1100             fjy0             = _mm_setzero_pd();
1101             fjz0             = _mm_setzero_pd();
1102
1103             /**************************
1104              * CALCULATE INTERACTIONS *
1105              **************************/
1106
1107             if (gmx_mm_any_lt(rsq00,rcutoff2))
1108             {
1109
1110             r00              = _mm_mul_pd(rsq00,rinv00);
1111
1112             /* Compute parameters for interactions between i and j atoms */
1113             gmx_mm_load_1pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,&c6_00,&c12_00);
1114
1115             /* LENNARD-JONES DISPERSION/REPULSION */
1116
1117             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
1118             vvdw6            = _mm_mul_pd(c6_00,rinvsix);
1119             vvdw12           = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
1120             vvdw             = _mm_sub_pd( _mm_mul_pd(vvdw12,one_twelfth) , _mm_mul_pd(vvdw6,one_sixth) );
1121             fvdw             = _mm_mul_pd(_mm_sub_pd(vvdw12,vvdw6),rinvsq00);
1122
1123             d                = _mm_sub_pd(r00,rswitch);
1124             d                = _mm_max_pd(d,_mm_setzero_pd());
1125             d2               = _mm_mul_pd(d,d);
1126             sw               = _mm_add_pd(one,_mm_mul_pd(d2,_mm_mul_pd(d,_mm_add_pd(swV3,_mm_mul_pd(d,_mm_add_pd(swV4,_mm_mul_pd(d,swV5)))))));
1127
1128             dsw              = _mm_mul_pd(d2,_mm_add_pd(swF2,_mm_mul_pd(d,_mm_add_pd(swF3,_mm_mul_pd(d,swF4)))));
1129
1130             /* Evaluate switch function */
1131             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1132             fvdw             = _mm_sub_pd( _mm_mul_pd(fvdw,sw) , _mm_mul_pd(rinv00,_mm_mul_pd(vvdw,dsw)) );
1133             cutoff_mask      = _mm_cmplt_pd(rsq00,rcutoff2);
1134
1135             fscal            = fvdw;
1136
1137             fscal            = _mm_and_pd(fscal,cutoff_mask);
1138
1139             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1140
1141             /* Calculate temporary vectorial force */
1142             tx               = _mm_mul_pd(fscal,dx00);
1143             ty               = _mm_mul_pd(fscal,dy00);
1144             tz               = _mm_mul_pd(fscal,dz00);
1145
1146             /* Update vectorial force */
1147             fix0             = _mm_add_pd(fix0,tx);
1148             fiy0             = _mm_add_pd(fiy0,ty);
1149             fiz0             = _mm_add_pd(fiz0,tz);
1150
1151             fjx0             = _mm_add_pd(fjx0,tx);
1152             fjy0             = _mm_add_pd(fjy0,ty);
1153             fjz0             = _mm_add_pd(fjz0,tz);
1154
1155             }
1156
1157             /**************************
1158              * CALCULATE INTERACTIONS *
1159              **************************/
1160
1161             if (gmx_mm_any_lt(rsq10,rcutoff2))
1162             {
1163
1164             /* Compute parameters for interactions between i and j atoms */
1165             qq10             = _mm_mul_pd(iq1,jq0);
1166
1167             /* REACTION-FIELD ELECTROSTATICS */
1168             felec            = _mm_mul_pd(qq10,_mm_sub_pd(_mm_mul_pd(rinv10,rinvsq10),krf2));
1169
1170             cutoff_mask      = _mm_cmplt_pd(rsq10,rcutoff2);
1171
1172             fscal            = felec;
1173
1174             fscal            = _mm_and_pd(fscal,cutoff_mask);
1175
1176             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1177
1178             /* Calculate temporary vectorial force */
1179             tx               = _mm_mul_pd(fscal,dx10);
1180             ty               = _mm_mul_pd(fscal,dy10);
1181             tz               = _mm_mul_pd(fscal,dz10);
1182
1183             /* Update vectorial force */
1184             fix1             = _mm_add_pd(fix1,tx);
1185             fiy1             = _mm_add_pd(fiy1,ty);
1186             fiz1             = _mm_add_pd(fiz1,tz);
1187
1188             fjx0             = _mm_add_pd(fjx0,tx);
1189             fjy0             = _mm_add_pd(fjy0,ty);
1190             fjz0             = _mm_add_pd(fjz0,tz);
1191
1192             }
1193
1194             /**************************
1195              * CALCULATE INTERACTIONS *
1196              **************************/
1197
1198             if (gmx_mm_any_lt(rsq20,rcutoff2))
1199             {
1200
1201             /* Compute parameters for interactions between i and j atoms */
1202             qq20             = _mm_mul_pd(iq2,jq0);
1203
1204             /* REACTION-FIELD ELECTROSTATICS */
1205             felec            = _mm_mul_pd(qq20,_mm_sub_pd(_mm_mul_pd(rinv20,rinvsq20),krf2));
1206
1207             cutoff_mask      = _mm_cmplt_pd(rsq20,rcutoff2);
1208
1209             fscal            = felec;
1210
1211             fscal            = _mm_and_pd(fscal,cutoff_mask);
1212
1213             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1214
1215             /* Calculate temporary vectorial force */
1216             tx               = _mm_mul_pd(fscal,dx20);
1217             ty               = _mm_mul_pd(fscal,dy20);
1218             tz               = _mm_mul_pd(fscal,dz20);
1219
1220             /* Update vectorial force */
1221             fix2             = _mm_add_pd(fix2,tx);
1222             fiy2             = _mm_add_pd(fiy2,ty);
1223             fiz2             = _mm_add_pd(fiz2,tz);
1224
1225             fjx0             = _mm_add_pd(fjx0,tx);
1226             fjy0             = _mm_add_pd(fjy0,ty);
1227             fjz0             = _mm_add_pd(fjz0,tz);
1228
1229             }
1230
1231             /**************************
1232              * CALCULATE INTERACTIONS *
1233              **************************/
1234
1235             if (gmx_mm_any_lt(rsq30,rcutoff2))
1236             {
1237
1238             /* Compute parameters for interactions between i and j atoms */
1239             qq30             = _mm_mul_pd(iq3,jq0);
1240
1241             /* REACTION-FIELD ELECTROSTATICS */
1242             felec            = _mm_mul_pd(qq30,_mm_sub_pd(_mm_mul_pd(rinv30,rinvsq30),krf2));
1243
1244             cutoff_mask      = _mm_cmplt_pd(rsq30,rcutoff2);
1245
1246             fscal            = felec;
1247
1248             fscal            = _mm_and_pd(fscal,cutoff_mask);
1249
1250             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1251
1252             /* Calculate temporary vectorial force */
1253             tx               = _mm_mul_pd(fscal,dx30);
1254             ty               = _mm_mul_pd(fscal,dy30);
1255             tz               = _mm_mul_pd(fscal,dz30);
1256
1257             /* Update vectorial force */
1258             fix3             = _mm_add_pd(fix3,tx);
1259             fiy3             = _mm_add_pd(fiy3,ty);
1260             fiz3             = _mm_add_pd(fiz3,tz);
1261
1262             fjx0             = _mm_add_pd(fjx0,tx);
1263             fjy0             = _mm_add_pd(fjy0,ty);
1264             fjz0             = _mm_add_pd(fjz0,tz);
1265
1266             }
1267
1268             gmx_mm_decrement_1rvec_1ptr_swizzle_pd(f+j_coord_offsetA,fjx0,fjy0,fjz0);
1269
1270             /* Inner loop uses 149 flops */
1271         }
1272
1273         /* End of innermost loop */
1274
1275         gmx_mm_update_iforce_4atom_swizzle_pd(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
1276                                               f+i_coord_offset,fshift+i_shift_offset);
1277
1278         /* Increment number of inner iterations */
1279         inneriter                  += j_index_end - j_index_start;
1280
1281         /* Outer loop uses 24 flops */
1282     }
1283
1284     /* Increment number of outer iterations */
1285     outeriter        += nri;
1286
1287     /* Update outer/inner flops */
1288
1289     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4_F,outeriter*24 + inneriter*149);
1290 }