src/gmxlib/nonbonded/nb_kernel_c/nb_kernel_ElecEwSh_VdwNone_GeomW4P1_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_VdwNone_GeomW4P1_VF_c
  35  * Electrostatics interaction: Ewald
  36  * VdW interaction:            None
  37  * Geometry:                   Water4-Particle
  38  * Calculate force/pot:        PotentialAndForce
  39  */
  40 void
  41 nb_kernel_ElecEwSh_VdwNone_GeomW4P1_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              vdwioffset1;
  57     real             ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
  58     int              vdwioffset2;
  59     real             ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
  60     int              vdwioffset3;
  61     real             ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
  62     int              vdwjidx0;
  63     real             jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
  64     real             dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10,cexp1_10,cexp2_10;
  65     real             dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20,cexp1_20,cexp2_20;
  66     real             dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30,cexp1_30,cexp2_30;
  67     real             velec,felec,velecsum,facel,crf,krf,krf2;
  68     real             *charge;
  69     int              ewitab;
  70     real             ewtabscale,eweps,sh_ewald,ewrt,ewtabhalfspace;
  71     real             *ewtab;
  72
  73     x                = xx[0];
  74     f                = ff[0];
  75
  76     nri              = nlist->nri;
  77     iinr             = nlist->iinr;
  78     jindex           = nlist->jindex;
  79     jjnr             = nlist->jjnr;
  80     shiftidx         = nlist->shift;
  81     gid              = nlist->gid;
  82     shiftvec         = fr->shift_vec[0];
  83     fshift           = fr->fshift[0];
  84     facel            = fr->epsfac;
  85     charge           = mdatoms->chargeA;
  86
  87     sh_ewald         = fr->ic->sh_ewald;
  88     ewtab            = fr->ic->tabq_coul_FDV0;
  89     ewtabscale       = fr->ic->tabq_scale;
  90     ewtabhalfspace   = 0.5/ewtabscale;
  91
  92     /* Setup water-specific parameters */
  93     inr              = nlist->iinr[0];
  94     iq1              = facel*charge[inr+1];
  95     iq2              = facel*charge[inr+2];
  96     iq3              = facel*charge[inr+3];
  97
  98     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
  99     rcutoff          = fr->rcoulomb;
 100     rcutoff2         = rcutoff*rcutoff;
 101
 102     outeriter        = 0;
 103     inneriter        = 0;
 104
 105     /* Start outer loop over neighborlists */
 106     for(iidx=0; iidx<nri; iidx++)
 107     {
 108         /* Load shift vector for this list */
 109         i_shift_offset   = DIM*shiftidx[iidx];
 110         shX              = shiftvec[i_shift_offset+XX];
 111         shY              = shiftvec[i_shift_offset+YY];
 112         shZ              = shiftvec[i_shift_offset+ZZ];
 113
 114         /* Load limits for loop over neighbors */
 115         j_index_start    = jindex[iidx];
 116         j_index_end      = jindex[iidx+1];
 117
 118         /* Get outer coordinate index */
 119         inr              = iinr[iidx];
 120         i_coord_offset   = DIM*inr;
 121
 122         /* Load i particle coords and add shift vector */
 123         ix1              = shX + x[i_coord_offset+DIM*1+XX];
 124         iy1              = shY + x[i_coord_offset+DIM*1+YY];
 125         iz1              = shZ + x[i_coord_offset+DIM*1+ZZ];
 126         ix2              = shX + x[i_coord_offset+DIM*2+XX];
 127         iy2              = shY + x[i_coord_offset+DIM*2+YY];
 128         iz2              = shZ + x[i_coord_offset+DIM*2+ZZ];
 129         ix3              = shX + x[i_coord_offset+DIM*3+XX];
 130         iy3              = shY + x[i_coord_offset+DIM*3+YY];
 131         iz3              = shZ + x[i_coord_offset+DIM*3+ZZ];
 132
 133         fix1             = 0.0;
 134         fiy1             = 0.0;
 135         fiz1             = 0.0;
 136         fix2             = 0.0;
 137         fiy2             = 0.0;
 138         fiz2             = 0.0;
 139         fix3             = 0.0;
 140         fiy3             = 0.0;
 141         fiz3             = 0.0;
 142
 143         /* Reset potential sums */
 144         velecsum         = 0.0;
 145
 146         /* Start inner kernel loop */
 147         for(jidx=j_index_start; jidx<j_index_end; jidx++)
 148         {
 149             /* Get j neighbor index, and coordinate index */
 150             jnr              = jjnr[jidx];
 151             j_coord_offset   = DIM*jnr;
 152
 153             /* load j atom coordinates */
 154             jx0              = x[j_coord_offset+DIM*0+XX];
 155             jy0              = x[j_coord_offset+DIM*0+YY];
 156             jz0              = x[j_coord_offset+DIM*0+ZZ];
 157
 158             /* Calculate displacement vector */
 159             dx10             = ix1 - jx0;
 160             dy10             = iy1 - jy0;
 161             dz10             = iz1 - jz0;
 162             dx20             = ix2 - jx0;
 163             dy20             = iy2 - jy0;
 164             dz20             = iz2 - jz0;
 165             dx30             = ix3 - jx0;
 166             dy30             = iy3 - jy0;
 167             dz30             = iz3 - jz0;
 168
 169             /* Calculate squared distance and things based on it */
 170             rsq10            = dx10*dx10+dy10*dy10+dz10*dz10;
 171             rsq20            = dx20*dx20+dy20*dy20+dz20*dz20;
 172             rsq30            = dx30*dx30+dy30*dy30+dz30*dz30;
 173
 174             rinv10           = gmx_invsqrt(rsq10);
 175             rinv20           = gmx_invsqrt(rsq20);
 176             rinv30           = gmx_invsqrt(rsq30);
 177
 178             rinvsq10         = rinv10*rinv10;
 179             rinvsq20         = rinv20*rinv20;
 180             rinvsq30         = rinv30*rinv30;
 181
 182             /* Load parameters for j particles */
 183             jq0              = charge[jnr+0];
 184
 185             /**************************
 186              * CALCULATE INTERACTIONS *
 187              **************************/
 188
 189             if (rsq10<rcutoff2)
 190             {
 191
 192             r10              = rsq10*rinv10;
 193
 194             qq10             = iq1*jq0;
 195
 196             /* EWALD ELECTROSTATICS */
 197
 198             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 199             ewrt             = r10*ewtabscale;
 200             ewitab           = ewrt;
 201             eweps            = ewrt-ewitab;
 202             ewitab           = 4*ewitab;
 203             felec            = ewtab[ewitab]+eweps*ewtab[ewitab+1];
 204             velec            = qq10*((rinv10-sh_ewald)-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
 205             felec            = qq10*rinv10*(rinvsq10-felec);
 206
 207             /* Update potential sums from outer loop */
 208             velecsum        += velec;
 209
 210             fscal            = felec;
 211
 212             /* Calculate temporary vectorial force */
 213             tx               = fscal*dx10;
 214             ty               = fscal*dy10;
 215             tz               = fscal*dz10;
 216
 217             /* Update vectorial force */
 218             fix1            += tx;
 219             fiy1            += ty;
 220             fiz1            += tz;
 221             f[j_coord_offset+DIM*0+XX] -= tx;
 222             f[j_coord_offset+DIM*0+YY] -= ty;
 223             f[j_coord_offset+DIM*0+ZZ] -= tz;
 224
 225             }
 226
 227             /**************************
 228              * CALCULATE INTERACTIONS *
 229              **************************/
 230
 231             if (rsq20<rcutoff2)
 232             {
 233
 234             r20              = rsq20*rinv20;
 235
 236             qq20             = iq2*jq0;
 237
 238             /* EWALD ELECTROSTATICS */
 239
 240             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 241             ewrt             = r20*ewtabscale;
 242             ewitab           = ewrt;
 243             eweps            = ewrt-ewitab;
 244             ewitab           = 4*ewitab;
 245             felec            = ewtab[ewitab]+eweps*ewtab[ewitab+1];
 246             velec            = qq20*((rinv20-sh_ewald)-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
 247             felec            = qq20*rinv20*(rinvsq20-felec);
 248
 249             /* Update potential sums from outer loop */
 250             velecsum        += velec;
 251
 252             fscal            = felec;
 253
 254             /* Calculate temporary vectorial force */
 255             tx               = fscal*dx20;
 256             ty               = fscal*dy20;
 257             tz               = fscal*dz20;
 258
 259             /* Update vectorial force */
 260             fix2            += tx;
 261             fiy2            += ty;
 262             fiz2            += tz;
 263             f[j_coord_offset+DIM*0+XX] -= tx;
 264             f[j_coord_offset+DIM*0+YY] -= ty;
 265             f[j_coord_offset+DIM*0+ZZ] -= tz;
 266
 267             }
 268
 269             /**************************
 270              * CALCULATE INTERACTIONS *
 271              **************************/
 272
 273             if (rsq30<rcutoff2)
 274             {
 275
 276             r30              = rsq30*rinv30;
 277
 278             qq30             = iq3*jq0;
 279
 280             /* EWALD ELECTROSTATICS */
 281
 282             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 283             ewrt             = r30*ewtabscale;
 284             ewitab           = ewrt;
 285             eweps            = ewrt-ewitab;
 286             ewitab           = 4*ewitab;
 287             felec            = ewtab[ewitab]+eweps*ewtab[ewitab+1];
 288             velec            = qq30*((rinv30-sh_ewald)-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
 289             felec            = qq30*rinv30*(rinvsq30-felec);
 290
 291             /* Update potential sums from outer loop */
 292             velecsum        += velec;
 293
 294             fscal            = felec;
 295
 296             /* Calculate temporary vectorial force */
 297             tx               = fscal*dx30;
 298             ty               = fscal*dy30;
 299             tz               = fscal*dz30;
 300
 301             /* Update vectorial force */
 302             fix3            += tx;
 303             fiy3            += ty;
 304             fiz3            += tz;
 305             f[j_coord_offset+DIM*0+XX] -= tx;
 306             f[j_coord_offset+DIM*0+YY] -= ty;
 307             f[j_coord_offset+DIM*0+ZZ] -= tz;
 308
 309             }
 310
 311             /* Inner loop uses 126 flops */
 312         }
 313         /* End of innermost loop */
 314
 315         tx = ty = tz = 0;
 316         f[i_coord_offset+DIM*1+XX] += fix1;
 317         f[i_coord_offset+DIM*1+YY] += fiy1;
 318         f[i_coord_offset+DIM*1+ZZ] += fiz1;
 319         tx                         += fix1;
 320         ty                         += fiy1;
 321         tz                         += fiz1;
 322         f[i_coord_offset+DIM*2+XX] += fix2;
 323         f[i_coord_offset+DIM*2+YY] += fiy2;
 324         f[i_coord_offset+DIM*2+ZZ] += fiz2;
 325         tx                         += fix2;
 326         ty                         += fiy2;
 327         tz                         += fiz2;
 328         f[i_coord_offset+DIM*3+XX] += fix3;
 329         f[i_coord_offset+DIM*3+YY] += fiy3;
 330         f[i_coord_offset+DIM*3+ZZ] += fiz3;
 331         tx                         += fix3;
 332         ty                         += fiy3;
 333         tz                         += fiz3;
 334         fshift[i_shift_offset+XX]  += tx;
 335         fshift[i_shift_offset+YY]  += ty;
 336         fshift[i_shift_offset+ZZ]  += tz;
 337
 338         ggid                        = gid[iidx];
 339         /* Update potential energies */
 340         kernel_data->energygrp_elec[ggid] += velecsum;
 341
 342         /* Increment number of inner iterations */
 343         inneriter                  += j_index_end - j_index_start;
 344
 345         /* Outer loop uses 31 flops */
 346     }
 347
 348     /* Increment number of outer iterations */
 349     outeriter        += nri;
 350
 351     /* Update outer/inner flops */
 352
 353     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_W4_VF,outeriter*31 + inneriter*126);
 354 }
 355 /*
 356  * Gromacs nonbonded kernel:   nb_kernel_ElecEwSh_VdwNone_GeomW4P1_F_c
 357  * Electrostatics interaction: Ewald
 358  * VdW interaction:            None
 359  * Geometry:                   Water4-Particle
 360  * Calculate force/pot:        Force
 361  */
 362 void
 363 nb_kernel_ElecEwSh_VdwNone_GeomW4P1_F_c
 364                     (t_nblist * gmx_restrict                nlist,
 365                      rvec * gmx_restrict                    xx,
 366                      rvec * gmx_restrict                    ff,
 367                      t_forcerec * gmx_restrict              fr,
 368                      t_mdatoms * gmx_restrict               mdatoms,
 369                      nb_kernel_data_t * gmx_restrict        kernel_data,
 370                      t_nrnb * gmx_restrict                  nrnb)
 371 {
 372     int              i_shift_offset,i_coord_offset,j_coord_offset;
 373     int              j_index_start,j_index_end;
 374     int              nri,inr,ggid,iidx,jidx,jnr,outeriter,inneriter;
 375     real             shX,shY,shZ,tx,ty,tz,fscal,rcutoff,rcutoff2;
 376     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
 377     real             *shiftvec,*fshift,*x,*f;
 378     int              vdwioffset1;
 379     real             ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
 380     int              vdwioffset2;
 381     real             ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
 382     int              vdwioffset3;
 383     real             ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
 384     int              vdwjidx0;
 385     real             jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
 386     real             dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10,cexp1_10,cexp2_10;
 387     real             dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20,cexp1_20,cexp2_20;
 388     real             dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30,cexp1_30,cexp2_30;
 389     real             velec,felec,velecsum,facel,crf,krf,krf2;
 390     real             *charge;
 391     int              ewitab;
 392     real             ewtabscale,eweps,sh_ewald,ewrt,ewtabhalfspace;
 393     real             *ewtab;
 394
 395     x                = xx[0];
 396     f                = ff[0];
 397
 398     nri              = nlist->nri;
 399     iinr             = nlist->iinr;
 400     jindex           = nlist->jindex;
 401     jjnr             = nlist->jjnr;
 402     shiftidx         = nlist->shift;
 403     gid              = nlist->gid;
 404     shiftvec         = fr->shift_vec[0];
 405     fshift           = fr->fshift[0];
 406     facel            = fr->epsfac;
 407     charge           = mdatoms->chargeA;
 408
 409     sh_ewald         = fr->ic->sh_ewald;
 410     ewtab            = fr->ic->tabq_coul_F;
 411     ewtabscale       = fr->ic->tabq_scale;
 412     ewtabhalfspace   = 0.5/ewtabscale;
 413
 414     /* Setup water-specific parameters */
 415     inr              = nlist->iinr[0];
 416     iq1              = facel*charge[inr+1];
 417     iq2              = facel*charge[inr+2];
 418     iq3              = facel*charge[inr+3];
 419
 420     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
 421     rcutoff          = fr->rcoulomb;
 422     rcutoff2         = rcutoff*rcutoff;
 423
 424     outeriter        = 0;
 425     inneriter        = 0;
 426
 427     /* Start outer loop over neighborlists */
 428     for(iidx=0; iidx<nri; iidx++)
 429     {
 430         /* Load shift vector for this list */
 431         i_shift_offset   = DIM*shiftidx[iidx];
 432         shX              = shiftvec[i_shift_offset+XX];
 433         shY              = shiftvec[i_shift_offset+YY];
 434         shZ              = shiftvec[i_shift_offset+ZZ];
 435
 436         /* Load limits for loop over neighbors */
 437         j_index_start    = jindex[iidx];
 438         j_index_end      = jindex[iidx+1];
 439
 440         /* Get outer coordinate index */
 441         inr              = iinr[iidx];
 442         i_coord_offset   = DIM*inr;
 443
 444         /* Load i particle coords and add shift vector */
 445         ix1              = shX + x[i_coord_offset+DIM*1+XX];
 446         iy1              = shY + x[i_coord_offset+DIM*1+YY];
 447         iz1              = shZ + x[i_coord_offset+DIM*1+ZZ];
 448         ix2              = shX + x[i_coord_offset+DIM*2+XX];
 449         iy2              = shY + x[i_coord_offset+DIM*2+YY];
 450         iz2              = shZ + x[i_coord_offset+DIM*2+ZZ];
 451         ix3              = shX + x[i_coord_offset+DIM*3+XX];
 452         iy3              = shY + x[i_coord_offset+DIM*3+YY];
 453         iz3              = shZ + x[i_coord_offset+DIM*3+ZZ];
 454
 455         fix1             = 0.0;
 456         fiy1             = 0.0;
 457         fiz1             = 0.0;
 458         fix2             = 0.0;
 459         fiy2             = 0.0;
 460         fiz2             = 0.0;
 461         fix3             = 0.0;
 462         fiy3             = 0.0;
 463         fiz3             = 0.0;
 464
 465         /* Start inner kernel loop */
 466         for(jidx=j_index_start; jidx<j_index_end; jidx++)
 467         {
 468             /* Get j neighbor index, and coordinate index */
 469             jnr              = jjnr[jidx];
 470             j_coord_offset   = DIM*jnr;
 471
 472             /* load j atom coordinates */
 473             jx0              = x[j_coord_offset+DIM*0+XX];
 474             jy0              = x[j_coord_offset+DIM*0+YY];
 475             jz0              = x[j_coord_offset+DIM*0+ZZ];
 476
 477             /* Calculate displacement vector */
 478             dx10             = ix1 - jx0;
 479             dy10             = iy1 - jy0;
 480             dz10             = iz1 - jz0;
 481             dx20             = ix2 - jx0;
 482             dy20             = iy2 - jy0;
 483             dz20             = iz2 - jz0;
 484             dx30             = ix3 - jx0;
 485             dy30             = iy3 - jy0;
 486             dz30             = iz3 - jz0;
 487
 488             /* Calculate squared distance and things based on it */
 489             rsq10            = dx10*dx10+dy10*dy10+dz10*dz10;
 490             rsq20            = dx20*dx20+dy20*dy20+dz20*dz20;
 491             rsq30            = dx30*dx30+dy30*dy30+dz30*dz30;
 492
 493             rinv10           = gmx_invsqrt(rsq10);
 494             rinv20           = gmx_invsqrt(rsq20);
 495             rinv30           = gmx_invsqrt(rsq30);
 496
 497             rinvsq10         = rinv10*rinv10;
 498             rinvsq20         = rinv20*rinv20;
 499             rinvsq30         = rinv30*rinv30;
 500
 501             /* Load parameters for j particles */
 502             jq0              = charge[jnr+0];
 503
 504             /**************************
 505              * CALCULATE INTERACTIONS *
 506              **************************/
 507
 508             if (rsq10<rcutoff2)
 509             {
 510
 511             r10              = rsq10*rinv10;
 512
 513             qq10             = iq1*jq0;
 514
 515             /* EWALD ELECTROSTATICS */
 516
 517             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 518             ewrt             = r10*ewtabscale;
 519             ewitab           = ewrt;
 520             eweps            = ewrt-ewitab;
 521             felec            = (1.0-eweps)*ewtab[ewitab]+eweps*ewtab[ewitab+1];
 522             felec            = qq10*rinv10*(rinvsq10-felec);
 523
 524             fscal            = felec;
 525
 526             /* Calculate temporary vectorial force */
 527             tx               = fscal*dx10;
 528             ty               = fscal*dy10;
 529             tz               = fscal*dz10;
 530
 531             /* Update vectorial force */
 532             fix1            += tx;
 533             fiy1            += ty;
 534             fiz1            += tz;
 535             f[j_coord_offset+DIM*0+XX] -= tx;
 536             f[j_coord_offset+DIM*0+YY] -= ty;
 537             f[j_coord_offset+DIM*0+ZZ] -= tz;
 538
 539             }
 540
 541             /**************************
 542              * CALCULATE INTERACTIONS *
 543              **************************/
 544
 545             if (rsq20<rcutoff2)
 546             {
 547
 548             r20              = rsq20*rinv20;
 549
 550             qq20             = iq2*jq0;
 551
 552             /* EWALD ELECTROSTATICS */
 553
 554             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 555             ewrt             = r20*ewtabscale;
 556             ewitab           = ewrt;
 557             eweps            = ewrt-ewitab;
 558             felec            = (1.0-eweps)*ewtab[ewitab]+eweps*ewtab[ewitab+1];
 559             felec            = qq20*rinv20*(rinvsq20-felec);
 560
 561             fscal            = felec;
 562
 563             /* Calculate temporary vectorial force */
 564             tx               = fscal*dx20;
 565             ty               = fscal*dy20;
 566             tz               = fscal*dz20;
 567
 568             /* Update vectorial force */
 569             fix2            += tx;
 570             fiy2            += ty;
 571             fiz2            += tz;
 572             f[j_coord_offset+DIM*0+XX] -= tx;
 573             f[j_coord_offset+DIM*0+YY] -= ty;
 574             f[j_coord_offset+DIM*0+ZZ] -= tz;
 575
 576             }
 577
 578             /**************************
 579              * CALCULATE INTERACTIONS *
 580              **************************/
 581
 582             if (rsq30<rcutoff2)
 583             {
 584
 585             r30              = rsq30*rinv30;
 586
 587             qq30             = iq3*jq0;
 588
 589             /* EWALD ELECTROSTATICS */
 590
 591             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 592             ewrt             = r30*ewtabscale;
 593             ewitab           = ewrt;
 594             eweps            = ewrt-ewitab;
 595             felec            = (1.0-eweps)*ewtab[ewitab]+eweps*ewtab[ewitab+1];
 596             felec            = qq30*rinv30*(rinvsq30-felec);
 597
 598             fscal            = felec;
 599
 600             /* Calculate temporary vectorial force */
 601             tx               = fscal*dx30;
 602             ty               = fscal*dy30;
 603             tz               = fscal*dz30;
 604
 605             /* Update vectorial force */
 606             fix3            += tx;
 607             fiy3            += ty;
 608             fiz3            += tz;
 609             f[j_coord_offset+DIM*0+XX] -= tx;
 610             f[j_coord_offset+DIM*0+YY] -= ty;
 611             f[j_coord_offset+DIM*0+ZZ] -= tz;
 612
 613             }
 614
 615             /* Inner loop uses 102 flops */
 616         }
 617         /* End of innermost loop */
 618
 619         tx = ty = tz = 0;
 620         f[i_coord_offset+DIM*1+XX] += fix1;
 621         f[i_coord_offset+DIM*1+YY] += fiy1;
 622         f[i_coord_offset+DIM*1+ZZ] += fiz1;
 623         tx                         += fix1;
 624         ty                         += fiy1;
 625         tz                         += fiz1;
 626         f[i_coord_offset+DIM*2+XX] += fix2;
 627         f[i_coord_offset+DIM*2+YY] += fiy2;
 628         f[i_coord_offset+DIM*2+ZZ] += fiz2;
 629         tx                         += fix2;
 630         ty                         += fiy2;
 631         tz                         += fiz2;
 632         f[i_coord_offset+DIM*3+XX] += fix3;
 633         f[i_coord_offset+DIM*3+YY] += fiy3;
 634         f[i_coord_offset+DIM*3+ZZ] += fiz3;
 635         tx                         += fix3;
 636         ty                         += fiy3;
 637         tz                         += fiz3;
 638         fshift[i_shift_offset+XX]  += tx;
 639         fshift[i_shift_offset+YY]  += ty;
 640         fshift[i_shift_offset+ZZ]  += tz;
 641
 642         /* Increment number of inner iterations */
 643         inneriter                  += j_index_end - j_index_start;
 644
 645         /* Outer loop uses 30 flops */
 646     }
 647
 648     /* Increment number of outer iterations */
 649     outeriter        += nri;
 650
 651     /* Update outer/inner flops */
 652
 653     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_W4_F,outeriter*30 + inneriter*102);
 654 }