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

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