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