src/gmxlib/nonbonded/nb_kernel_c/nb_kernel_ElecEwSh_VdwLJSh_GeomW4W4_c.c

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