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