src/gmxlib/nonbonded/nb_kernel_c/nb_kernel_ElecEw_VdwBham_GeomW4W4_c.c

   1 /*
   2  * Note: this file was generated by the Gromacs c 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 /*
  34  * Gromacs nonbonded kernel:   nb_kernel_ElecEw_VdwBham_GeomW4W4_VF_c
  35  * Electrostatics interaction: Ewald
  36  * VdW interaction:            Buckingham
  37  * Geometry:                   Water4-Water4
  38  * Calculate force/pot:        PotentialAndForce
  39  */
  40 void
  41 nb_kernel_ElecEw_VdwBham_GeomW4W4_VF_c
  42                     (t_nblist * gmx_restrict                nlist,
  43                      rvec * gmx_restrict                    xx,
  44                      rvec * gmx_restrict                    ff,
  45                      t_forcerec * gmx_restrict              fr,
  46                      t_mdatoms * gmx_restrict               mdatoms,
  47                      nb_kernel_data_t * gmx_restrict        kernel_data,
  48                      t_nrnb * gmx_restrict                  nrnb)
  49 {
  50     int              i_shift_offset,i_coord_offset,j_coord_offset;
  51     int              j_index_start,j_index_end;
  52     int              nri,inr,ggid,iidx,jidx,jnr,outeriter,inneriter;
  53     real             shX,shY,shZ,tx,ty,tz,fscal,rcutoff,rcutoff2;
  54     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
  55     real             *shiftvec,*fshift,*x,*f;
  56     int              vdwioffset0;
  57     real             ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
  58     int              vdwioffset1;
  59     real             ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
  60     int              vdwioffset2;
  61     real             ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
  62     int              vdwioffset3;
  63     real             ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
  64     int              vdwjidx0;
  65     real             jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
  66     int              vdwjidx1;
  67     real             jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
  68     int              vdwjidx2;
  69     real             jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
  70     int              vdwjidx3;
  71     real             jx3,jy3,jz3,fjx3,fjy3,fjz3,jq3,isaj3;
  72     real             dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00,cexp1_00,cexp2_00;
  73     real             dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11,cexp1_11,cexp2_11;
  74     real             dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12,cexp1_12,cexp2_12;
  75     real             dx13,dy13,dz13,rsq13,rinv13,rinvsq13,r13,qq13,c6_13,c12_13,cexp1_13,cexp2_13;
  76     real             dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21,cexp1_21,cexp2_21;
  77     real             dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22,cexp1_22,cexp2_22;
  78     real             dx23,dy23,dz23,rsq23,rinv23,rinvsq23,r23,qq23,c6_23,c12_23,cexp1_23,cexp2_23;
  79     real             dx31,dy31,dz31,rsq31,rinv31,rinvsq31,r31,qq31,c6_31,c12_31,cexp1_31,cexp2_31;
  80     real             dx32,dy32,dz32,rsq32,rinv32,rinvsq32,r32,qq32,c6_32,c12_32,cexp1_32,cexp2_32;
  81     real             dx33,dy33,dz33,rsq33,rinv33,rinvsq33,r33,qq33,c6_33,c12_33,cexp1_33,cexp2_33;
  82     real             velec,felec,velecsum,facel,crf,krf,krf2;
  83     real             *charge;
  84     int              nvdwtype;
  85     real             rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,br,vvdwexp,sh_vdw_invrcut6;
  86     int              *vdwtype;
  87     real             *vdwparam;
  88     int              ewitab;
  89     real             ewtabscale,eweps,sh_ewald,ewrt,ewtabhalfspace;
  90     real             *ewtab;
  91
  92     x                = xx[0];
  93     f                = ff[0];
  94
  95     nri              = nlist->nri;
  96     iinr             = nlist->iinr;
  97     jindex           = nlist->jindex;
  98     jjnr             = nlist->jjnr;
  99     shiftidx         = nlist->shift;
 100     gid              = nlist->gid;
 101     shiftvec         = fr->shift_vec[0];
 102     fshift           = fr->fshift[0];
 103     facel            = fr->epsfac;
 104     charge           = mdatoms->chargeA;
 105     nvdwtype         = fr->ntype;
 106     vdwparam         = fr->nbfp;
 107     vdwtype          = mdatoms->typeA;
 108
 109     sh_ewald         = fr->ic->sh_ewald;
 110     ewtab            = fr->ic->tabq_coul_FDV0;
 111     ewtabscale       = fr->ic->tabq_scale;
 112     ewtabhalfspace   = 0.5/ewtabscale;
 113
 114     /* Setup water-specific parameters */
 115     inr              = nlist->iinr[0];
 116     iq1              = facel*charge[inr+1];
 117     iq2              = facel*charge[inr+2];
 118     iq3              = facel*charge[inr+3];
 119     vdwioffset0      = 3*nvdwtype*vdwtype[inr+0];
 120
 121     jq1              = charge[inr+1];
 122     jq2              = charge[inr+2];
 123     jq3              = charge[inr+3];
 124     vdwjidx0         = 3*vdwtype[inr+0];
 125     c6_00            = vdwparam[vdwioffset0+vdwjidx0];
 126     cexp1_00         = vdwparam[vdwioffset0+vdwjidx0+1];
 127     cexp2_00         = vdwparam[vdwioffset0+vdwjidx0+2];
 128     qq11             = iq1*jq1;
 129     qq12             = iq1*jq2;
 130     qq13             = iq1*jq3;
 131     qq21             = iq2*jq1;
 132     qq22             = iq2*jq2;
 133     qq23             = iq2*jq3;
 134     qq31             = iq3*jq1;
 135     qq32             = iq3*jq2;
 136     qq33             = iq3*jq3;
 137
 138     outeriter        = 0;
 139     inneriter        = 0;
 140
 141     /* Start outer loop over neighborlists */
 142     for(iidx=0; iidx<nri; iidx++)
 143     {
 144         /* Load shift vector for this list */
 145         i_shift_offset   = DIM*shiftidx[iidx];
 146         shX              = shiftvec[i_shift_offset+XX];
 147         shY              = shiftvec[i_shift_offset+YY];
 148         shZ              = shiftvec[i_shift_offset+ZZ];
 149
 150         /* Load limits for loop over neighbors */
 151         j_index_start    = jindex[iidx];
 152         j_index_end      = jindex[iidx+1];
 153
 154         /* Get outer coordinate index */
 155         inr              = iinr[iidx];
 156         i_coord_offset   = DIM*inr;
 157
 158         /* Load i particle coords and add shift vector */
 159         ix0              = shX + x[i_coord_offset+DIM*0+XX];
 160         iy0              = shY + x[i_coord_offset+DIM*0+YY];
 161         iz0              = shZ + x[i_coord_offset+DIM*0+ZZ];
 162         ix1              = shX + x[i_coord_offset+DIM*1+XX];
 163         iy1              = shY + x[i_coord_offset+DIM*1+YY];
 164         iz1              = shZ + x[i_coord_offset+DIM*1+ZZ];
 165         ix2              = shX + x[i_coord_offset+DIM*2+XX];
 166         iy2              = shY + x[i_coord_offset+DIM*2+YY];
 167         iz2              = shZ + x[i_coord_offset+DIM*2+ZZ];
 168         ix3              = shX + x[i_coord_offset+DIM*3+XX];
 169         iy3              = shY + x[i_coord_offset+DIM*3+YY];
 170         iz3              = shZ + x[i_coord_offset+DIM*3+ZZ];
 171
 172         fix0             = 0.0;
 173         fiy0             = 0.0;
 174         fiz0             = 0.0;
 175         fix1             = 0.0;
 176         fiy1             = 0.0;
 177         fiz1             = 0.0;
 178         fix2             = 0.0;
 179         fiy2             = 0.0;
 180         fiz2             = 0.0;
 181         fix3             = 0.0;
 182         fiy3             = 0.0;
 183         fiz3             = 0.0;
 184
 185         /* Reset potential sums */
 186         velecsum         = 0.0;
 187         vvdwsum          = 0.0;
 188
 189         /* Start inner kernel loop */
 190         for(jidx=j_index_start; jidx<j_index_end; jidx++)
 191         {
 192             /* Get j neighbor index, and coordinate index */
 193             jnr              = jjnr[jidx];
 194             j_coord_offset   = DIM*jnr;
 195
 196             /* load j atom coordinates */
 197             jx0              = x[j_coord_offset+DIM*0+XX];
 198             jy0              = x[j_coord_offset+DIM*0+YY];
 199             jz0              = x[j_coord_offset+DIM*0+ZZ];
 200             jx1              = x[j_coord_offset+DIM*1+XX];
 201             jy1              = x[j_coord_offset+DIM*1+YY];
 202             jz1              = x[j_coord_offset+DIM*1+ZZ];
 203             jx2              = x[j_coord_offset+DIM*2+XX];
 204             jy2              = x[j_coord_offset+DIM*2+YY];
 205             jz2              = x[j_coord_offset+DIM*2+ZZ];
 206             jx3              = x[j_coord_offset+DIM*3+XX];
 207             jy3              = x[j_coord_offset+DIM*3+YY];
 208             jz3              = x[j_coord_offset+DIM*3+ZZ];
 209
 210             /* Calculate displacement vector */
 211             dx00             = ix0 - jx0;
 212             dy00             = iy0 - jy0;
 213             dz00             = iz0 - jz0;
 214             dx11             = ix1 - jx1;
 215             dy11             = iy1 - jy1;
 216             dz11             = iz1 - jz1;
 217             dx12             = ix1 - jx2;
 218             dy12             = iy1 - jy2;
 219             dz12             = iz1 - jz2;
 220             dx13             = ix1 - jx3;
 221             dy13             = iy1 - jy3;
 222             dz13             = iz1 - jz3;
 223             dx21             = ix2 - jx1;
 224             dy21             = iy2 - jy1;
 225             dz21             = iz2 - jz1;
 226             dx22             = ix2 - jx2;
 227             dy22             = iy2 - jy2;
 228             dz22             = iz2 - jz2;
 229             dx23             = ix2 - jx3;
 230             dy23             = iy2 - jy3;
 231             dz23             = iz2 - jz3;
 232             dx31             = ix3 - jx1;
 233             dy31             = iy3 - jy1;
 234             dz31             = iz3 - jz1;
 235             dx32             = ix3 - jx2;
 236             dy32             = iy3 - jy2;
 237             dz32             = iz3 - jz2;
 238             dx33             = ix3 - jx3;
 239             dy33             = iy3 - jy3;
 240             dz33             = iz3 - jz3;
 241
 242             /* Calculate squared distance and things based on it */
 243             rsq00            = dx00*dx00+dy00*dy00+dz00*dz00;
 244             rsq11            = dx11*dx11+dy11*dy11+dz11*dz11;
 245             rsq12            = dx12*dx12+dy12*dy12+dz12*dz12;
 246             rsq13            = dx13*dx13+dy13*dy13+dz13*dz13;
 247             rsq21            = dx21*dx21+dy21*dy21+dz21*dz21;
 248             rsq22            = dx22*dx22+dy22*dy22+dz22*dz22;
 249             rsq23            = dx23*dx23+dy23*dy23+dz23*dz23;
 250             rsq31            = dx31*dx31+dy31*dy31+dz31*dz31;
 251             rsq32            = dx32*dx32+dy32*dy32+dz32*dz32;
 252             rsq33            = dx33*dx33+dy33*dy33+dz33*dz33;
 253
 254             rinv00           = gmx_invsqrt(rsq00);
 255             rinv11           = gmx_invsqrt(rsq11);
 256             rinv12           = gmx_invsqrt(rsq12);
 257             rinv13           = gmx_invsqrt(rsq13);
 258             rinv21           = gmx_invsqrt(rsq21);
 259             rinv22           = gmx_invsqrt(rsq22);
 260             rinv23           = gmx_invsqrt(rsq23);
 261             rinv31           = gmx_invsqrt(rsq31);
 262             rinv32           = gmx_invsqrt(rsq32);
 263             rinv33           = gmx_invsqrt(rsq33);
 264
 265             rinvsq00         = rinv00*rinv00;
 266             rinvsq11         = rinv11*rinv11;
 267             rinvsq12         = rinv12*rinv12;
 268             rinvsq13         = rinv13*rinv13;
 269             rinvsq21         = rinv21*rinv21;
 270             rinvsq22         = rinv22*rinv22;
 271             rinvsq23         = rinv23*rinv23;
 272             rinvsq31         = rinv31*rinv31;
 273             rinvsq32         = rinv32*rinv32;
 274             rinvsq33         = rinv33*rinv33;
 275
 276             /**************************
 277              * CALCULATE INTERACTIONS *
 278              **************************/
 279
 280             r00              = rsq00*rinv00;
 281
 282             /* BUCKINGHAM DISPERSION/REPULSION */
 283             rinvsix          = rinvsq00*rinvsq00*rinvsq00;
 284             vvdw6            = c6_00*rinvsix;
 285             br               = cexp2_00*r00;
 286             vvdwexp          = cexp1_00*exp(-br);
 287             vvdw             = vvdwexp - vvdw6*(1.0/6.0);
 288             fvdw             = (br*vvdwexp-vvdw6)*rinvsq00;
 289
 290             /* Update potential sums from outer loop */
 291             vvdwsum         += vvdw;
 292
 293             fscal            = fvdw;
 294
 295             /* Calculate temporary vectorial force */
 296             tx               = fscal*dx00;
 297             ty               = fscal*dy00;
 298             tz               = fscal*dz00;
 299
 300             /* Update vectorial force */
 301             fix0            += tx;
 302             fiy0            += ty;
 303             fiz0            += tz;
 304             f[j_coord_offset+DIM*0+XX] -= tx;
 305             f[j_coord_offset+DIM*0+YY] -= ty;
 306             f[j_coord_offset+DIM*0+ZZ] -= tz;
 307
 308             /**************************
 309              * CALCULATE INTERACTIONS *
 310              **************************/
 311
 312             r11              = rsq11*rinv11;
 313
 314             /* EWALD ELECTROSTATICS */
 315
 316             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 317             ewrt             = r11*ewtabscale;
 318             ewitab           = ewrt;
 319             eweps            = ewrt-ewitab;
 320             ewitab           = 4*ewitab;
 321             felec            = ewtab[ewitab]+eweps*ewtab[ewitab+1];
 322             velec            = qq11*(rinv11-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
 323             felec            = qq11*rinv11*(rinvsq11-felec);
 324
 325             /* Update potential sums from outer loop */
 326             velecsum        += velec;
 327
 328             fscal            = felec;
 329
 330             /* Calculate temporary vectorial force */
 331             tx               = fscal*dx11;
 332             ty               = fscal*dy11;
 333             tz               = fscal*dz11;
 334
 335             /* Update vectorial force */
 336             fix1            += tx;
 337             fiy1            += ty;
 338             fiz1            += tz;
 339             f[j_coord_offset+DIM*1+XX] -= tx;
 340             f[j_coord_offset+DIM*1+YY] -= ty;
 341             f[j_coord_offset+DIM*1+ZZ] -= tz;
 342
 343             /**************************
 344              * CALCULATE INTERACTIONS *
 345              **************************/
 346
 347             r12              = rsq12*rinv12;
 348
 349             /* EWALD ELECTROSTATICS */
 350
 351             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 352             ewrt             = r12*ewtabscale;
 353             ewitab           = ewrt;
 354             eweps            = ewrt-ewitab;
 355             ewitab           = 4*ewitab;
 356             felec            = ewtab[ewitab]+eweps*ewtab[ewitab+1];
 357             velec            = qq12*(rinv12-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
 358             felec            = qq12*rinv12*(rinvsq12-felec);
 359
 360             /* Update potential sums from outer loop */
 361             velecsum        += velec;
 362
 363             fscal            = felec;
 364
 365             /* Calculate temporary vectorial force */
 366             tx               = fscal*dx12;
 367             ty               = fscal*dy12;
 368             tz               = fscal*dz12;
 369
 370             /* Update vectorial force */
 371             fix1            += tx;
 372             fiy1            += ty;
 373             fiz1            += tz;
 374             f[j_coord_offset+DIM*2+XX] -= tx;
 375             f[j_coord_offset+DIM*2+YY] -= ty;
 376             f[j_coord_offset+DIM*2+ZZ] -= tz;
 377
 378             /**************************
 379              * CALCULATE INTERACTIONS *
 380              **************************/
 381
 382             r13              = rsq13*rinv13;
 383
 384             /* EWALD ELECTROSTATICS */
 385
 386             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 387             ewrt             = r13*ewtabscale;
 388             ewitab           = ewrt;
 389             eweps            = ewrt-ewitab;
 390             ewitab           = 4*ewitab;
 391             felec            = ewtab[ewitab]+eweps*ewtab[ewitab+1];
 392             velec            = qq13*(rinv13-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
 393             felec            = qq13*rinv13*(rinvsq13-felec);
 394
 395             /* Update potential sums from outer loop */
 396             velecsum        += velec;
 397
 398             fscal            = felec;
 399
 400             /* Calculate temporary vectorial force */
 401             tx               = fscal*dx13;
 402             ty               = fscal*dy13;
 403             tz               = fscal*dz13;
 404
 405             /* Update vectorial force */
 406             fix1            += tx;
 407             fiy1            += ty;
 408             fiz1            += tz;
 409             f[j_coord_offset+DIM*3+XX] -= tx;
 410             f[j_coord_offset+DIM*3+YY] -= ty;
 411             f[j_coord_offset+DIM*3+ZZ] -= tz;
 412
 413             /**************************
 414              * CALCULATE INTERACTIONS *
 415              **************************/
 416
 417             r21              = rsq21*rinv21;
 418
 419             /* EWALD ELECTROSTATICS */
 420
 421             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 422             ewrt             = r21*ewtabscale;
 423             ewitab           = ewrt;
 424             eweps            = ewrt-ewitab;
 425             ewitab           = 4*ewitab;
 426             felec            = ewtab[ewitab]+eweps*ewtab[ewitab+1];
 427             velec            = qq21*(rinv21-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
 428             felec            = qq21*rinv21*(rinvsq21-felec);
 429
 430             /* Update potential sums from outer loop */
 431             velecsum        += velec;
 432
 433             fscal            = felec;
 434
 435             /* Calculate temporary vectorial force */
 436             tx               = fscal*dx21;
 437             ty               = fscal*dy21;
 438             tz               = fscal*dz21;
 439
 440             /* Update vectorial force */
 441             fix2            += tx;
 442             fiy2            += ty;
 443             fiz2            += tz;
 444             f[j_coord_offset+DIM*1+XX] -= tx;
 445             f[j_coord_offset+DIM*1+YY] -= ty;
 446             f[j_coord_offset+DIM*1+ZZ] -= tz;
 447
 448             /**************************
 449              * CALCULATE INTERACTIONS *
 450              **************************/
 451
 452             r22              = rsq22*rinv22;
 453
 454             /* EWALD ELECTROSTATICS */
 455
 456             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 457             ewrt             = r22*ewtabscale;
 458             ewitab           = ewrt;
 459             eweps            = ewrt-ewitab;
 460             ewitab           = 4*ewitab;
 461             felec            = ewtab[ewitab]+eweps*ewtab[ewitab+1];
 462             velec            = qq22*(rinv22-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
 463             felec            = qq22*rinv22*(rinvsq22-felec);
 464
 465             /* Update potential sums from outer loop */
 466             velecsum        += velec;
 467
 468             fscal            = felec;
 469
 470             /* Calculate temporary vectorial force */
 471             tx               = fscal*dx22;
 472             ty               = fscal*dy22;
 473             tz               = fscal*dz22;
 474
 475             /* Update vectorial force */
 476             fix2            += tx;
 477             fiy2            += ty;
 478             fiz2            += tz;
 479             f[j_coord_offset+DIM*2+XX] -= tx;
 480             f[j_coord_offset+DIM*2+YY] -= ty;
 481             f[j_coord_offset+DIM*2+ZZ] -= tz;
 482
 483             /**************************
 484              * CALCULATE INTERACTIONS *
 485              **************************/
 486
 487             r23              = rsq23*rinv23;
 488
 489             /* EWALD ELECTROSTATICS */
 490
 491             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 492             ewrt             = r23*ewtabscale;
 493             ewitab           = ewrt;
 494             eweps            = ewrt-ewitab;
 495             ewitab           = 4*ewitab;
 496             felec            = ewtab[ewitab]+eweps*ewtab[ewitab+1];
 497             velec            = qq23*(rinv23-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
 498             felec            = qq23*rinv23*(rinvsq23-felec);
 499
 500             /* Update potential sums from outer loop */
 501             velecsum        += velec;
 502
 503             fscal            = felec;
 504
 505             /* Calculate temporary vectorial force */
 506             tx               = fscal*dx23;
 507             ty               = fscal*dy23;
 508             tz               = fscal*dz23;
 509
 510             /* Update vectorial force */
 511             fix2            += tx;
 512             fiy2            += ty;
 513             fiz2            += tz;
 514             f[j_coord_offset+DIM*3+XX] -= tx;
 515             f[j_coord_offset+DIM*3+YY] -= ty;
 516             f[j_coord_offset+DIM*3+ZZ] -= tz;
 517
 518             /**************************
 519              * CALCULATE INTERACTIONS *
 520              **************************/
 521
 522             r31              = rsq31*rinv31;
 523
 524             /* EWALD ELECTROSTATICS */
 525
 526             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 527             ewrt             = r31*ewtabscale;
 528             ewitab           = ewrt;
 529             eweps            = ewrt-ewitab;
 530             ewitab           = 4*ewitab;
 531             felec            = ewtab[ewitab]+eweps*ewtab[ewitab+1];
 532             velec            = qq31*(rinv31-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
 533             felec            = qq31*rinv31*(rinvsq31-felec);
 534
 535             /* Update potential sums from outer loop */
 536             velecsum        += velec;
 537
 538             fscal            = felec;
 539
 540             /* Calculate temporary vectorial force */
 541             tx               = fscal*dx31;
 542             ty               = fscal*dy31;
 543             tz               = fscal*dz31;
 544
 545             /* Update vectorial force */
 546             fix3            += tx;
 547             fiy3            += ty;
 548             fiz3            += tz;
 549             f[j_coord_offset+DIM*1+XX] -= tx;
 550             f[j_coord_offset+DIM*1+YY] -= ty;
 551             f[j_coord_offset+DIM*1+ZZ] -= tz;
 552
 553             /**************************
 554              * CALCULATE INTERACTIONS *
 555              **************************/
 556
 557             r32              = rsq32*rinv32;
 558
 559             /* EWALD ELECTROSTATICS */
 560
 561             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 562             ewrt             = r32*ewtabscale;
 563             ewitab           = ewrt;
 564             eweps            = ewrt-ewitab;
 565             ewitab           = 4*ewitab;
 566             felec            = ewtab[ewitab]+eweps*ewtab[ewitab+1];
 567             velec            = qq32*(rinv32-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
 568             felec            = qq32*rinv32*(rinvsq32-felec);
 569
 570             /* Update potential sums from outer loop */
 571             velecsum        += velec;
 572
 573             fscal            = felec;
 574
 575             /* Calculate temporary vectorial force */
 576             tx               = fscal*dx32;
 577             ty               = fscal*dy32;
 578             tz               = fscal*dz32;
 579
 580             /* Update vectorial force */
 581             fix3            += tx;
 582             fiy3            += ty;
 583             fiz3            += tz;
 584             f[j_coord_offset+DIM*2+XX] -= tx;
 585             f[j_coord_offset+DIM*2+YY] -= ty;
 586             f[j_coord_offset+DIM*2+ZZ] -= tz;
 587
 588             /**************************
 589              * CALCULATE INTERACTIONS *
 590              **************************/
 591
 592             r33              = rsq33*rinv33;
 593
 594             /* EWALD ELECTROSTATICS */
 595
 596             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 597             ewrt             = r33*ewtabscale;
 598             ewitab           = ewrt;
 599             eweps            = ewrt-ewitab;
 600             ewitab           = 4*ewitab;
 601             felec            = ewtab[ewitab]+eweps*ewtab[ewitab+1];
 602             velec            = qq33*(rinv33-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
 603             felec            = qq33*rinv33*(rinvsq33-felec);
 604
 605             /* Update potential sums from outer loop */
 606             velecsum        += velec;
 607
 608             fscal            = felec;
 609
 610             /* Calculate temporary vectorial force */
 611             tx               = fscal*dx33;
 612             ty               = fscal*dy33;
 613             tz               = fscal*dz33;
 614
 615             /* Update vectorial force */
 616             fix3            += tx;
 617             fiy3            += ty;
 618             fiz3            += tz;
 619             f[j_coord_offset+DIM*3+XX] -= tx;
 620             f[j_coord_offset+DIM*3+YY] -= ty;
 621             f[j_coord_offset+DIM*3+ZZ] -= tz;
 622
 623             /* Inner loop uses 421 flops */
 624         }
 625         /* End of innermost loop */
 626
 627         tx = ty = tz = 0;
 628         f[i_coord_offset+DIM*0+XX] += fix0;
 629         f[i_coord_offset+DIM*0+YY] += fiy0;
 630         f[i_coord_offset+DIM*0+ZZ] += fiz0;
 631         tx                         += fix0;
 632         ty                         += fiy0;
 633         tz                         += fiz0;
 634         f[i_coord_offset+DIM*1+XX] += fix1;
 635         f[i_coord_offset+DIM*1+YY] += fiy1;
 636         f[i_coord_offset+DIM*1+ZZ] += fiz1;
 637         tx                         += fix1;
 638         ty                         += fiy1;
 639         tz                         += fiz1;
 640         f[i_coord_offset+DIM*2+XX] += fix2;
 641         f[i_coord_offset+DIM*2+YY] += fiy2;
 642         f[i_coord_offset+DIM*2+ZZ] += fiz2;
 643         tx                         += fix2;
 644         ty                         += fiy2;
 645         tz                         += fiz2;
 646         f[i_coord_offset+DIM*3+XX] += fix3;
 647         f[i_coord_offset+DIM*3+YY] += fiy3;
 648         f[i_coord_offset+DIM*3+ZZ] += fiz3;
 649         tx                         += fix3;
 650         ty                         += fiy3;
 651         tz                         += fiz3;
 652         fshift[i_shift_offset+XX]  += tx;
 653         fshift[i_shift_offset+YY]  += ty;
 654         fshift[i_shift_offset+ZZ]  += tz;
 655
 656         ggid                        = gid[iidx];
 657         /* Update potential energies */
 658         kernel_data->energygrp_elec[ggid] += velecsum;
 659         kernel_data->energygrp_vdw[ggid] += vvdwsum;
 660
 661         /* Increment number of inner iterations */
 662         inneriter                  += j_index_end - j_index_start;
 663
 664         /* Outer loop uses 41 flops */
 665     }
 666
 667     /* Increment number of outer iterations */
 668     outeriter        += nri;
 669
 670     /* Update outer/inner flops */
 671
 672     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4W4_VF,outeriter*41 + inneriter*421);
 673 }
 674 /*
 675  * Gromacs nonbonded kernel:   nb_kernel_ElecEw_VdwBham_GeomW4W4_F_c
 676  * Electrostatics interaction: Ewald
 677  * VdW interaction:            Buckingham
 678  * Geometry:                   Water4-Water4
 679  * Calculate force/pot:        Force
 680  */
 681 void
 682 nb_kernel_ElecEw_VdwBham_GeomW4W4_F_c
 683                     (t_nblist * gmx_restrict                nlist,
 684                      rvec * gmx_restrict                    xx,
 685                      rvec * gmx_restrict                    ff,
 686                      t_forcerec * gmx_restrict              fr,
 687                      t_mdatoms * gmx_restrict               mdatoms,
 688                      nb_kernel_data_t * gmx_restrict        kernel_data,
 689                      t_nrnb * gmx_restrict                  nrnb)
 690 {
 691     int              i_shift_offset,i_coord_offset,j_coord_offset;
 692     int              j_index_start,j_index_end;
 693     int              nri,inr,ggid,iidx,jidx,jnr,outeriter,inneriter;
 694     real             shX,shY,shZ,tx,ty,tz,fscal,rcutoff,rcutoff2;
 695     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
 696     real             *shiftvec,*fshift,*x,*f;
 697     int              vdwioffset0;
 698     real             ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
 699     int              vdwioffset1;
 700     real             ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
 701     int              vdwioffset2;
 702     real             ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
 703     int              vdwioffset3;
 704     real             ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
 705     int              vdwjidx0;
 706     real             jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
 707     int              vdwjidx1;
 708     real             jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
 709     int              vdwjidx2;
 710     real             jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
 711     int              vdwjidx3;
 712     real             jx3,jy3,jz3,fjx3,fjy3,fjz3,jq3,isaj3;
 713     real             dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00,cexp1_00,cexp2_00;
 714     real             dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11,cexp1_11,cexp2_11;
 715     real             dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12,cexp1_12,cexp2_12;
 716     real             dx13,dy13,dz13,rsq13,rinv13,rinvsq13,r13,qq13,c6_13,c12_13,cexp1_13,cexp2_13;
 717     real             dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21,cexp1_21,cexp2_21;
 718     real             dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22,cexp1_22,cexp2_22;
 719     real             dx23,dy23,dz23,rsq23,rinv23,rinvsq23,r23,qq23,c6_23,c12_23,cexp1_23,cexp2_23;
 720     real             dx31,dy31,dz31,rsq31,rinv31,rinvsq31,r31,qq31,c6_31,c12_31,cexp1_31,cexp2_31;
 721     real             dx32,dy32,dz32,rsq32,rinv32,rinvsq32,r32,qq32,c6_32,c12_32,cexp1_32,cexp2_32;
 722     real             dx33,dy33,dz33,rsq33,rinv33,rinvsq33,r33,qq33,c6_33,c12_33,cexp1_33,cexp2_33;
 723     real             velec,felec,velecsum,facel,crf,krf,krf2;
 724     real             *charge;
 725     int              nvdwtype;
 726     real             rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,br,vvdwexp,sh_vdw_invrcut6;
 727     int              *vdwtype;
 728     real             *vdwparam;
 729     int              ewitab;
 730     real             ewtabscale,eweps,sh_ewald,ewrt,ewtabhalfspace;
 731     real             *ewtab;
 732
 733     x                = xx[0];
 734     f                = ff[0];
 735
 736     nri              = nlist->nri;
 737     iinr             = nlist->iinr;
 738     jindex           = nlist->jindex;
 739     jjnr             = nlist->jjnr;
 740     shiftidx         = nlist->shift;
 741     gid              = nlist->gid;
 742     shiftvec         = fr->shift_vec[0];
 743     fshift           = fr->fshift[0];
 744     facel            = fr->epsfac;
 745     charge           = mdatoms->chargeA;
 746     nvdwtype         = fr->ntype;
 747     vdwparam         = fr->nbfp;
 748     vdwtype          = mdatoms->typeA;
 749
 750     sh_ewald         = fr->ic->sh_ewald;
 751     ewtab            = fr->ic->tabq_coul_F;
 752     ewtabscale       = fr->ic->tabq_scale;
 753     ewtabhalfspace   = 0.5/ewtabscale;
 754
 755     /* Setup water-specific parameters */
 756     inr              = nlist->iinr[0];
 757     iq1              = facel*charge[inr+1];
 758     iq2              = facel*charge[inr+2];
 759     iq3              = facel*charge[inr+3];
 760     vdwioffset0      = 3*nvdwtype*vdwtype[inr+0];
 761
 762     jq1              = charge[inr+1];
 763     jq2              = charge[inr+2];
 764     jq3              = charge[inr+3];
 765     vdwjidx0         = 3*vdwtype[inr+0];
 766     c6_00            = vdwparam[vdwioffset0+vdwjidx0];
 767     cexp1_00         = vdwparam[vdwioffset0+vdwjidx0+1];
 768     cexp2_00         = vdwparam[vdwioffset0+vdwjidx0+2];
 769     qq11             = iq1*jq1;
 770     qq12             = iq1*jq2;
 771     qq13             = iq1*jq3;
 772     qq21             = iq2*jq1;
 773     qq22             = iq2*jq2;
 774     qq23             = iq2*jq3;
 775     qq31             = iq3*jq1;
 776     qq32             = iq3*jq2;
 777     qq33             = iq3*jq3;
 778
 779     outeriter        = 0;
 780     inneriter        = 0;
 781
 782     /* Start outer loop over neighborlists */
 783     for(iidx=0; iidx<nri; iidx++)
 784     {
 785         /* Load shift vector for this list */
 786         i_shift_offset   = DIM*shiftidx[iidx];
 787         shX              = shiftvec[i_shift_offset+XX];
 788         shY              = shiftvec[i_shift_offset+YY];
 789         shZ              = shiftvec[i_shift_offset+ZZ];
 790
 791         /* Load limits for loop over neighbors */
 792         j_index_start    = jindex[iidx];
 793         j_index_end      = jindex[iidx+1];
 794
 795         /* Get outer coordinate index */
 796         inr              = iinr[iidx];
 797         i_coord_offset   = DIM*inr;
 798
 799         /* Load i particle coords and add shift vector */
 800         ix0              = shX + x[i_coord_offset+DIM*0+XX];
 801         iy0              = shY + x[i_coord_offset+DIM*0+YY];
 802         iz0              = shZ + x[i_coord_offset+DIM*0+ZZ];
 803         ix1              = shX + x[i_coord_offset+DIM*1+XX];
 804         iy1              = shY + x[i_coord_offset+DIM*1+YY];
 805         iz1              = shZ + x[i_coord_offset+DIM*1+ZZ];
 806         ix2              = shX + x[i_coord_offset+DIM*2+XX];
 807         iy2              = shY + x[i_coord_offset+DIM*2+YY];
 808         iz2              = shZ + x[i_coord_offset+DIM*2+ZZ];
 809         ix3              = shX + x[i_coord_offset+DIM*3+XX];
 810         iy3              = shY + x[i_coord_offset+DIM*3+YY];
 811         iz3              = shZ + x[i_coord_offset+DIM*3+ZZ];
 812
 813         fix0             = 0.0;
 814         fiy0             = 0.0;
 815         fiz0             = 0.0;
 816         fix1             = 0.0;
 817         fiy1             = 0.0;
 818         fiz1             = 0.0;
 819         fix2             = 0.0;
 820         fiy2             = 0.0;
 821         fiz2             = 0.0;
 822         fix3             = 0.0;
 823         fiy3             = 0.0;
 824         fiz3             = 0.0;
 825
 826         /* Start inner kernel loop */
 827         for(jidx=j_index_start; jidx<j_index_end; jidx++)
 828         {
 829             /* Get j neighbor index, and coordinate index */
 830             jnr              = jjnr[jidx];
 831             j_coord_offset   = DIM*jnr;
 832
 833             /* load j atom coordinates */
 834             jx0              = x[j_coord_offset+DIM*0+XX];
 835             jy0              = x[j_coord_offset+DIM*0+YY];
 836             jz0              = x[j_coord_offset+DIM*0+ZZ];
 837             jx1              = x[j_coord_offset+DIM*1+XX];
 838             jy1              = x[j_coord_offset+DIM*1+YY];
 839             jz1              = x[j_coord_offset+DIM*1+ZZ];
 840             jx2              = x[j_coord_offset+DIM*2+XX];
 841             jy2              = x[j_coord_offset+DIM*2+YY];
 842             jz2              = x[j_coord_offset+DIM*2+ZZ];
 843             jx3              = x[j_coord_offset+DIM*3+XX];
 844             jy3              = x[j_coord_offset+DIM*3+YY];
 845             jz3              = x[j_coord_offset+DIM*3+ZZ];
 846
 847             /* Calculate displacement vector */
 848             dx00             = ix0 - jx0;
 849             dy00             = iy0 - jy0;
 850             dz00             = iz0 - jz0;
 851             dx11             = ix1 - jx1;
 852             dy11             = iy1 - jy1;
 853             dz11             = iz1 - jz1;
 854             dx12             = ix1 - jx2;
 855             dy12             = iy1 - jy2;
 856             dz12             = iz1 - jz2;
 857             dx13             = ix1 - jx3;
 858             dy13             = iy1 - jy3;
 859             dz13             = iz1 - jz3;
 860             dx21             = ix2 - jx1;
 861             dy21             = iy2 - jy1;
 862             dz21             = iz2 - jz1;
 863             dx22             = ix2 - jx2;
 864             dy22             = iy2 - jy2;
 865             dz22             = iz2 - jz2;
 866             dx23             = ix2 - jx3;
 867             dy23             = iy2 - jy3;
 868             dz23             = iz2 - jz3;
 869             dx31             = ix3 - jx1;
 870             dy31             = iy3 - jy1;
 871             dz31             = iz3 - jz1;
 872             dx32             = ix3 - jx2;
 873             dy32             = iy3 - jy2;
 874             dz32             = iz3 - jz2;
 875             dx33             = ix3 - jx3;
 876             dy33             = iy3 - jy3;
 877             dz33             = iz3 - jz3;
 878
 879             /* Calculate squared distance and things based on it */
 880             rsq00            = dx00*dx00+dy00*dy00+dz00*dz00;
 881             rsq11            = dx11*dx11+dy11*dy11+dz11*dz11;
 882             rsq12            = dx12*dx12+dy12*dy12+dz12*dz12;
 883             rsq13            = dx13*dx13+dy13*dy13+dz13*dz13;
 884             rsq21            = dx21*dx21+dy21*dy21+dz21*dz21;
 885             rsq22            = dx22*dx22+dy22*dy22+dz22*dz22;
 886             rsq23            = dx23*dx23+dy23*dy23+dz23*dz23;
 887             rsq31            = dx31*dx31+dy31*dy31+dz31*dz31;
 888             rsq32            = dx32*dx32+dy32*dy32+dz32*dz32;
 889             rsq33            = dx33*dx33+dy33*dy33+dz33*dz33;
 890
 891             rinv00           = gmx_invsqrt(rsq00);
 892             rinv11           = gmx_invsqrt(rsq11);
 893             rinv12           = gmx_invsqrt(rsq12);
 894             rinv13           = gmx_invsqrt(rsq13);
 895             rinv21           = gmx_invsqrt(rsq21);
 896             rinv22           = gmx_invsqrt(rsq22);
 897             rinv23           = gmx_invsqrt(rsq23);
 898             rinv31           = gmx_invsqrt(rsq31);
 899             rinv32           = gmx_invsqrt(rsq32);
 900             rinv33           = gmx_invsqrt(rsq33);
 901
 902             rinvsq00         = rinv00*rinv00;
 903             rinvsq11         = rinv11*rinv11;
 904             rinvsq12         = rinv12*rinv12;
 905             rinvsq13         = rinv13*rinv13;
 906             rinvsq21         = rinv21*rinv21;
 907             rinvsq22         = rinv22*rinv22;
 908             rinvsq23         = rinv23*rinv23;
 909             rinvsq31         = rinv31*rinv31;
 910             rinvsq32         = rinv32*rinv32;
 911             rinvsq33         = rinv33*rinv33;
 912
 913             /**************************
 914              * CALCULATE INTERACTIONS *
 915              **************************/
 916
 917             r00              = rsq00*rinv00;
 918
 919             /* BUCKINGHAM DISPERSION/REPULSION */
 920             rinvsix          = rinvsq00*rinvsq00*rinvsq00;
 921             vvdw6            = c6_00*rinvsix;
 922             br               = cexp2_00*r00;
 923             vvdwexp          = cexp1_00*exp(-br);
 924             fvdw             = (br*vvdwexp-vvdw6)*rinvsq00;
 925
 926             fscal            = fvdw;
 927
 928             /* Calculate temporary vectorial force */
 929             tx               = fscal*dx00;
 930             ty               = fscal*dy00;
 931             tz               = fscal*dz00;
 932
 933             /* Update vectorial force */
 934             fix0            += tx;
 935             fiy0            += ty;
 936             fiz0            += tz;
 937             f[j_coord_offset+DIM*0+XX] -= tx;
 938             f[j_coord_offset+DIM*0+YY] -= ty;
 939             f[j_coord_offset+DIM*0+ZZ] -= tz;
 940
 941             /**************************
 942              * CALCULATE INTERACTIONS *
 943              **************************/
 944
 945             r11              = rsq11*rinv11;
 946
 947             /* EWALD ELECTROSTATICS */
 948
 949             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 950             ewrt             = r11*ewtabscale;
 951             ewitab           = ewrt;
 952             eweps            = ewrt-ewitab;
 953             felec            = (1.0-eweps)*ewtab[ewitab]+eweps*ewtab[ewitab+1];
 954             felec            = qq11*rinv11*(rinvsq11-felec);
 955
 956             fscal            = felec;
 957
 958             /* Calculate temporary vectorial force */
 959             tx               = fscal*dx11;
 960             ty               = fscal*dy11;
 961             tz               = fscal*dz11;
 962
 963             /* Update vectorial force */
 964             fix1            += tx;
 965             fiy1            += ty;
 966             fiz1            += tz;
 967             f[j_coord_offset+DIM*1+XX] -= tx;
 968             f[j_coord_offset+DIM*1+YY] -= ty;
 969             f[j_coord_offset+DIM*1+ZZ] -= tz;
 970
 971             /**************************
 972              * CALCULATE INTERACTIONS *
 973              **************************/
 974
 975             r12              = rsq12*rinv12;
 976
 977             /* EWALD ELECTROSTATICS */
 978
 979             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 980             ewrt             = r12*ewtabscale;
 981             ewitab           = ewrt;
 982             eweps            = ewrt-ewitab;
 983             felec            = (1.0-eweps)*ewtab[ewitab]+eweps*ewtab[ewitab+1];
 984             felec            = qq12*rinv12*(rinvsq12-felec);
 985
 986             fscal            = felec;
 987
 988             /* Calculate temporary vectorial force */
 989             tx               = fscal*dx12;
 990             ty               = fscal*dy12;
 991             tz               = fscal*dz12;
 992
 993             /* Update vectorial force */
 994             fix1            += tx;
 995             fiy1            += ty;
 996             fiz1            += tz;
 997             f[j_coord_offset+DIM*2+XX] -= tx;
 998             f[j_coord_offset+DIM*2+YY] -= ty;
 999             f[j_coord_offset+DIM*2+ZZ] -= tz;
1000
1001             /**************************
1002              * CALCULATE INTERACTIONS *
1003              **************************/
1004
1005             r13              = rsq13*rinv13;
1006
1007             /* EWALD ELECTROSTATICS */
1008
1009             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1010             ewrt             = r13*ewtabscale;
1011             ewitab           = ewrt;
1012             eweps            = ewrt-ewitab;
1013             felec            = (1.0-eweps)*ewtab[ewitab]+eweps*ewtab[ewitab+1];
1014             felec            = qq13*rinv13*(rinvsq13-felec);
1015
1016             fscal            = felec;
1017
1018             /* Calculate temporary vectorial force */
1019             tx               = fscal*dx13;
1020             ty               = fscal*dy13;
1021             tz               = fscal*dz13;
1022
1023             /* Update vectorial force */
1024             fix1            += tx;
1025             fiy1            += ty;
1026             fiz1            += tz;
1027             f[j_coord_offset+DIM*3+XX] -= tx;
1028             f[j_coord_offset+DIM*3+YY] -= ty;
1029             f[j_coord_offset+DIM*3+ZZ] -= tz;
1030
1031             /**************************
1032              * CALCULATE INTERACTIONS *
1033              **************************/
1034
1035             r21              = rsq21*rinv21;
1036
1037             /* EWALD ELECTROSTATICS */
1038
1039             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1040             ewrt             = r21*ewtabscale;
1041             ewitab           = ewrt;
1042             eweps            = ewrt-ewitab;
1043             felec            = (1.0-eweps)*ewtab[ewitab]+eweps*ewtab[ewitab+1];
1044             felec            = qq21*rinv21*(rinvsq21-felec);
1045
1046             fscal            = felec;
1047
1048             /* Calculate temporary vectorial force */
1049             tx               = fscal*dx21;
1050             ty               = fscal*dy21;
1051             tz               = fscal*dz21;
1052
1053             /* Update vectorial force */
1054             fix2            += tx;
1055             fiy2            += ty;
1056             fiz2            += tz;
1057             f[j_coord_offset+DIM*1+XX] -= tx;
1058             f[j_coord_offset+DIM*1+YY] -= ty;
1059             f[j_coord_offset+DIM*1+ZZ] -= tz;
1060
1061             /**************************
1062              * CALCULATE INTERACTIONS *
1063              **************************/
1064
1065             r22              = rsq22*rinv22;
1066
1067             /* EWALD ELECTROSTATICS */
1068
1069             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1070             ewrt             = r22*ewtabscale;
1071             ewitab           = ewrt;
1072             eweps            = ewrt-ewitab;
1073             felec            = (1.0-eweps)*ewtab[ewitab]+eweps*ewtab[ewitab+1];
1074             felec            = qq22*rinv22*(rinvsq22-felec);
1075
1076             fscal            = felec;
1077
1078             /* Calculate temporary vectorial force */
1079             tx               = fscal*dx22;
1080             ty               = fscal*dy22;
1081             tz               = fscal*dz22;
1082
1083             /* Update vectorial force */
1084             fix2            += tx;
1085             fiy2            += ty;
1086             fiz2            += tz;
1087             f[j_coord_offset+DIM*2+XX] -= tx;
1088             f[j_coord_offset+DIM*2+YY] -= ty;
1089             f[j_coord_offset+DIM*2+ZZ] -= tz;
1090
1091             /**************************
1092              * CALCULATE INTERACTIONS *
1093              **************************/
1094
1095             r23              = rsq23*rinv23;
1096
1097             /* EWALD ELECTROSTATICS */
1098
1099             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1100             ewrt             = r23*ewtabscale;
1101             ewitab           = ewrt;
1102             eweps            = ewrt-ewitab;
1103             felec            = (1.0-eweps)*ewtab[ewitab]+eweps*ewtab[ewitab+1];
1104             felec            = qq23*rinv23*(rinvsq23-felec);
1105
1106             fscal            = felec;
1107
1108             /* Calculate temporary vectorial force */
1109             tx               = fscal*dx23;
1110             ty               = fscal*dy23;
1111             tz               = fscal*dz23;
1112
1113             /* Update vectorial force */
1114             fix2            += tx;
1115             fiy2            += ty;
1116             fiz2            += tz;
1117             f[j_coord_offset+DIM*3+XX] -= tx;
1118             f[j_coord_offset+DIM*3+YY] -= ty;
1119             f[j_coord_offset+DIM*3+ZZ] -= tz;
1120
1121             /**************************
1122              * CALCULATE INTERACTIONS *
1123              **************************/
1124
1125             r31              = rsq31*rinv31;
1126
1127             /* EWALD ELECTROSTATICS */
1128
1129             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1130             ewrt             = r31*ewtabscale;
1131             ewitab           = ewrt;
1132             eweps            = ewrt-ewitab;
1133             felec            = (1.0-eweps)*ewtab[ewitab]+eweps*ewtab[ewitab+1];
1134             felec            = qq31*rinv31*(rinvsq31-felec);
1135
1136             fscal            = felec;
1137
1138             /* Calculate temporary vectorial force */
1139             tx               = fscal*dx31;
1140             ty               = fscal*dy31;
1141             tz               = fscal*dz31;
1142
1143             /* Update vectorial force */
1144             fix3            += tx;
1145             fiy3            += ty;
1146             fiz3            += tz;
1147             f[j_coord_offset+DIM*1+XX] -= tx;
1148             f[j_coord_offset+DIM*1+YY] -= ty;
1149             f[j_coord_offset+DIM*1+ZZ] -= tz;
1150
1151             /**************************
1152              * CALCULATE INTERACTIONS *
1153              **************************/
1154
1155             r32              = rsq32*rinv32;
1156
1157             /* EWALD ELECTROSTATICS */
1158
1159             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1160             ewrt             = r32*ewtabscale;
1161             ewitab           = ewrt;
1162             eweps            = ewrt-ewitab;
1163             felec            = (1.0-eweps)*ewtab[ewitab]+eweps*ewtab[ewitab+1];
1164             felec            = qq32*rinv32*(rinvsq32-felec);
1165
1166             fscal            = felec;
1167
1168             /* Calculate temporary vectorial force */
1169             tx               = fscal*dx32;
1170             ty               = fscal*dy32;
1171             tz               = fscal*dz32;
1172
1173             /* Update vectorial force */
1174             fix3            += tx;
1175             fiy3            += ty;
1176             fiz3            += tz;
1177             f[j_coord_offset+DIM*2+XX] -= tx;
1178             f[j_coord_offset+DIM*2+YY] -= ty;
1179             f[j_coord_offset+DIM*2+ZZ] -= tz;
1180
1181             /**************************
1182              * CALCULATE INTERACTIONS *
1183              **************************/
1184
1185             r33              = rsq33*rinv33;
1186
1187             /* EWALD ELECTROSTATICS */
1188
1189             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1190             ewrt             = r33*ewtabscale;
1191             ewitab           = ewrt;
1192             eweps            = ewrt-ewitab;
1193             felec            = (1.0-eweps)*ewtab[ewitab]+eweps*ewtab[ewitab+1];
1194             felec            = qq33*rinv33*(rinvsq33-felec);
1195
1196             fscal            = felec;
1197
1198             /* Calculate temporary vectorial force */
1199             tx               = fscal*dx33;
1200             ty               = fscal*dy33;
1201             tz               = fscal*dz33;
1202
1203             /* Update vectorial force */
1204             fix3            += tx;
1205             fiy3            += ty;
1206             fiz3            += tz;
1207             f[j_coord_offset+DIM*3+XX] -= tx;
1208             f[j_coord_offset+DIM*3+YY] -= ty;
1209             f[j_coord_offset+DIM*3+ZZ] -= tz;
1210
1211             /* Inner loop uses 355 flops */
1212         }
1213         /* End of innermost loop */
1214
1215         tx = ty = tz = 0;
1216         f[i_coord_offset+DIM*0+XX] += fix0;
1217         f[i_coord_offset+DIM*0+YY] += fiy0;
1218         f[i_coord_offset+DIM*0+ZZ] += fiz0;
1219         tx                         += fix0;
1220         ty                         += fiy0;
1221         tz                         += fiz0;
1222         f[i_coord_offset+DIM*1+XX] += fix1;
1223         f[i_coord_offset+DIM*1+YY] += fiy1;
1224         f[i_coord_offset+DIM*1+ZZ] += fiz1;
1225         tx                         += fix1;
1226         ty                         += fiy1;
1227         tz                         += fiz1;
1228         f[i_coord_offset+DIM*2+XX] += fix2;
1229         f[i_coord_offset+DIM*2+YY] += fiy2;
1230         f[i_coord_offset+DIM*2+ZZ] += fiz2;
1231         tx                         += fix2;
1232         ty                         += fiy2;
1233         tz                         += fiz2;
1234         f[i_coord_offset+DIM*3+XX] += fix3;
1235         f[i_coord_offset+DIM*3+YY] += fiy3;
1236         f[i_coord_offset+DIM*3+ZZ] += fiz3;
1237         tx                         += fix3;
1238         ty                         += fiy3;
1239         tz                         += fiz3;
1240         fshift[i_shift_offset+XX]  += tx;
1241         fshift[i_shift_offset+YY]  += ty;
1242         fshift[i_shift_offset+ZZ]  += tz;
1243
1244         /* Increment number of inner iterations */
1245         inneriter                  += j_index_end - j_index_start;
1246
1247         /* Outer loop uses 39 flops */
1248     }
1249
1250     /* Increment number of outer iterations */
1251     outeriter        += nri;
1252
1253     /* Update outer/inner flops */
1254
1255     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4W4_F,outeriter*39 + inneriter*355);
1256 }