src/gromacs/gmxlib/nonbonded/nb_kernel_sse2_double/nb_kernel_ElecRFCut_VdwNone_GeomW4W4_sse2_double.c

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