src/gromacs/gmxlib/nonbonded/nb_kernel_sse2_double/nb_kernel_ElecCoul_VdwLJ_GeomW3W3_sse2_double.cpp

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