src/gromacs/gmxlib/nonbonded/nb_kernel_avx_128_fma_double/nb_kernel_ElecRFCut_VdwLJSw_GeomW3W3_avx_128_fma_double.c

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