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