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