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