src/gromacs/gmxlib/nonbonded/nb_kernel_c/nb_kernel_ElecRFCut_VdwNone_GeomW3P1_c.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*
  36  * Note: this file was generated by the GROMACS c kernel generator.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "config.h"
  41
  42 #include <math.h>
  43
  44 #include "../nb_kernel.h"
  45 #include "gromacs/legacyheaders/types/simple.h"
  46 #include "gromacs/math/vec.h"
  47 #include "gromacs/legacyheaders/nrnb.h"
  48
  49 /*
  50  * Gromacs nonbonded kernel:   nb_kernel_ElecRFCut_VdwNone_GeomW3P1_VF_c
  51  * Electrostatics interaction: ReactionField
  52  * VdW interaction:            None
  53  * Geometry:                   Water3-Particle
  54  * Calculate force/pot:        PotentialAndForce
  55  */
  56 void
  57 nb_kernel_ElecRFCut_VdwNone_GeomW3P1_VF_c
  58                     (t_nblist                    * gmx_restrict       nlist,
  59                      rvec                        * gmx_restrict          xx,
  60                      rvec                        * gmx_restrict          ff,
  61                      t_forcerec                  * gmx_restrict          fr,
  62                      t_mdatoms                   * gmx_restrict     mdatoms,
  63                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
  64                      t_nrnb                      * gmx_restrict        nrnb)
  65 {
  66     int              i_shift_offset,i_coord_offset,j_coord_offset;
  67     int              j_index_start,j_index_end;
  68     int              nri,inr,ggid,iidx,jidx,jnr,outeriter,inneriter;
  69     real             shX,shY,shZ,tx,ty,tz,fscal,rcutoff,rcutoff2;
  70     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
  71     real             *shiftvec,*fshift,*x,*f;
  72     int              vdwioffset0;
  73     real             ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
  74     int              vdwioffset1;
  75     real             ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
  76     int              vdwioffset2;
  77     real             ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
  78     int              vdwjidx0;
  79     real             jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
  80     real             dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00,cexp1_00,cexp2_00;
  81     real             dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10,cexp1_10,cexp2_10;
  82     real             dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20,cexp1_20,cexp2_20;
  83     real             velec,felec,velecsum,facel,crf,krf,krf2;
  84     real             *charge;
  85
  86     x                = xx[0];
  87     f                = ff[0];
  88
  89     nri              = nlist->nri;
  90     iinr             = nlist->iinr;
  91     jindex           = nlist->jindex;
  92     jjnr             = nlist->jjnr;
  93     shiftidx         = nlist->shift;
  94     gid              = nlist->gid;
  95     shiftvec         = fr->shift_vec[0];
  96     fshift           = fr->fshift[0];
  97     facel            = fr->epsfac;
  98     charge           = mdatoms->chargeA;
  99     krf              = fr->ic->k_rf;
 100     krf2             = krf*2.0;
 101     crf              = fr->ic->c_rf;
 102
 103     /* Setup water-specific parameters */
 104     inr              = nlist->iinr[0];
 105     iq0              = facel*charge[inr+0];
 106     iq1              = facel*charge[inr+1];
 107     iq2              = facel*charge[inr+2];
 108
 109     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
 110     rcutoff          = fr->rcoulomb;
 111     rcutoff2         = rcutoff*rcutoff;
 112
 113     outeriter        = 0;
 114     inneriter        = 0;
 115
 116     /* Start outer loop over neighborlists */
 117     for(iidx=0; iidx<nri; iidx++)
 118     {
 119         /* Load shift vector for this list */
 120         i_shift_offset   = DIM*shiftidx[iidx];
 121         shX              = shiftvec[i_shift_offset+XX];
 122         shY              = shiftvec[i_shift_offset+YY];
 123         shZ              = shiftvec[i_shift_offset+ZZ];
 124
 125         /* Load limits for loop over neighbors */
 126         j_index_start    = jindex[iidx];
 127         j_index_end      = jindex[iidx+1];
 128
 129         /* Get outer coordinate index */
 130         inr              = iinr[iidx];
 131         i_coord_offset   = DIM*inr;
 132
 133         /* Load i particle coords and add shift vector */
 134         ix0              = shX + x[i_coord_offset+DIM*0+XX];
 135         iy0              = shY + x[i_coord_offset+DIM*0+YY];
 136         iz0              = shZ + x[i_coord_offset+DIM*0+ZZ];
 137         ix1              = shX + x[i_coord_offset+DIM*1+XX];
 138         iy1              = shY + x[i_coord_offset+DIM*1+YY];
 139         iz1              = shZ + x[i_coord_offset+DIM*1+ZZ];
 140         ix2              = shX + x[i_coord_offset+DIM*2+XX];
 141         iy2              = shY + x[i_coord_offset+DIM*2+YY];
 142         iz2              = shZ + x[i_coord_offset+DIM*2+ZZ];
 143
 144         fix0             = 0.0;
 145         fiy0             = 0.0;
 146         fiz0             = 0.0;
 147         fix1             = 0.0;
 148         fiy1             = 0.0;
 149         fiz1             = 0.0;
 150         fix2             = 0.0;
 151         fiy2             = 0.0;
 152         fiz2             = 0.0;
 153
 154         /* Reset potential sums */
 155         velecsum         = 0.0;
 156
 157         /* Start inner kernel loop */
 158         for(jidx=j_index_start; jidx<j_index_end; jidx++)
 159         {
 160             /* Get j neighbor index, and coordinate index */
 161             jnr              = jjnr[jidx];
 162             j_coord_offset   = DIM*jnr;
 163
 164             /* load j atom coordinates */
 165             jx0              = x[j_coord_offset+DIM*0+XX];
 166             jy0              = x[j_coord_offset+DIM*0+YY];
 167             jz0              = x[j_coord_offset+DIM*0+ZZ];
 168
 169             /* Calculate displacement vector */
 170             dx00             = ix0 - jx0;
 171             dy00             = iy0 - jy0;
 172             dz00             = iz0 - jz0;
 173             dx10             = ix1 - jx0;
 174             dy10             = iy1 - jy0;
 175             dz10             = iz1 - jz0;
 176             dx20             = ix2 - jx0;
 177             dy20             = iy2 - jy0;
 178             dz20             = iz2 - jz0;
 179
 180             /* Calculate squared distance and things based on it */
 181             rsq00            = dx00*dx00+dy00*dy00+dz00*dz00;
 182             rsq10            = dx10*dx10+dy10*dy10+dz10*dz10;
 183             rsq20            = dx20*dx20+dy20*dy20+dz20*dz20;
 184
 185             rinv00           = gmx_invsqrt(rsq00);
 186             rinv10           = gmx_invsqrt(rsq10);
 187             rinv20           = gmx_invsqrt(rsq20);
 188
 189             rinvsq00         = rinv00*rinv00;
 190             rinvsq10         = rinv10*rinv10;
 191             rinvsq20         = rinv20*rinv20;
 192
 193             /* Load parameters for j particles */
 194             jq0              = charge[jnr+0];
 195
 196             /**************************
 197              * CALCULATE INTERACTIONS *
 198              **************************/
 199
 200             if (rsq00<rcutoff2)
 201             {
 202
 203             qq00             = iq0*jq0;
 204
 205             /* REACTION-FIELD ELECTROSTATICS */
 206             velec            = qq00*(rinv00+krf*rsq00-crf);
 207             felec            = qq00*(rinv00*rinvsq00-krf2);
 208
 209             /* Update potential sums from outer loop */
 210             velecsum        += velec;
 211
 212             fscal            = felec;
 213
 214             /* Calculate temporary vectorial force */
 215             tx               = fscal*dx00;
 216             ty               = fscal*dy00;
 217             tz               = fscal*dz00;
 218
 219             /* Update vectorial force */
 220             fix0            += tx;
 221             fiy0            += ty;
 222             fiz0            += tz;
 223             f[j_coord_offset+DIM*0+XX] -= tx;
 224             f[j_coord_offset+DIM*0+YY] -= ty;
 225             f[j_coord_offset+DIM*0+ZZ] -= tz;
 226
 227             }
 228
 229             /**************************
 230              * CALCULATE INTERACTIONS *
 231              **************************/
 232
 233             if (rsq10<rcutoff2)
 234             {
 235
 236             qq10             = iq1*jq0;
 237
 238             /* REACTION-FIELD ELECTROSTATICS */
 239             velec            = qq10*(rinv10+krf*rsq10-crf);
 240             felec            = qq10*(rinv10*rinvsq10-krf2);
 241
 242             /* Update potential sums from outer loop */
 243             velecsum        += velec;
 244
 245             fscal            = felec;
 246
 247             /* Calculate temporary vectorial force */
 248             tx               = fscal*dx10;
 249             ty               = fscal*dy10;
 250             tz               = fscal*dz10;
 251
 252             /* Update vectorial force */
 253             fix1            += tx;
 254             fiy1            += ty;
 255             fiz1            += tz;
 256             f[j_coord_offset+DIM*0+XX] -= tx;
 257             f[j_coord_offset+DIM*0+YY] -= ty;
 258             f[j_coord_offset+DIM*0+ZZ] -= tz;
 259
 260             }
 261
 262             /**************************
 263              * CALCULATE INTERACTIONS *
 264              **************************/
 265
 266             if (rsq20<rcutoff2)
 267             {
 268
 269             qq20             = iq2*jq0;
 270
 271             /* REACTION-FIELD ELECTROSTATICS */
 272             velec            = qq20*(rinv20+krf*rsq20-crf);
 273             felec            = qq20*(rinv20*rinvsq20-krf2);
 274
 275             /* Update potential sums from outer loop */
 276             velecsum        += velec;
 277
 278             fscal            = felec;
 279
 280             /* Calculate temporary vectorial force */
 281             tx               = fscal*dx20;
 282             ty               = fscal*dy20;
 283             tz               = fscal*dz20;
 284
 285             /* Update vectorial force */
 286             fix2            += tx;
 287             fiy2            += ty;
 288             fiz2            += tz;
 289             f[j_coord_offset+DIM*0+XX] -= tx;
 290             f[j_coord_offset+DIM*0+YY] -= ty;
 291             f[j_coord_offset+DIM*0+ZZ] -= tz;
 292
 293             }
 294
 295             /* Inner loop uses 96 flops */
 296         }
 297         /* End of innermost loop */
 298
 299         tx = ty = tz = 0;
 300         f[i_coord_offset+DIM*0+XX] += fix0;
 301         f[i_coord_offset+DIM*0+YY] += fiy0;
 302         f[i_coord_offset+DIM*0+ZZ] += fiz0;
 303         tx                         += fix0;
 304         ty                         += fiy0;
 305         tz                         += fiz0;
 306         f[i_coord_offset+DIM*1+XX] += fix1;
 307         f[i_coord_offset+DIM*1+YY] += fiy1;
 308         f[i_coord_offset+DIM*1+ZZ] += fiz1;
 309         tx                         += fix1;
 310         ty                         += fiy1;
 311         tz                         += fiz1;
 312         f[i_coord_offset+DIM*2+XX] += fix2;
 313         f[i_coord_offset+DIM*2+YY] += fiy2;
 314         f[i_coord_offset+DIM*2+ZZ] += fiz2;
 315         tx                         += fix2;
 316         ty                         += fiy2;
 317         tz                         += fiz2;
 318         fshift[i_shift_offset+XX]  += tx;
 319         fshift[i_shift_offset+YY]  += ty;
 320         fshift[i_shift_offset+ZZ]  += tz;
 321
 322         ggid                        = gid[iidx];
 323         /* Update potential energies */
 324         kernel_data->energygrp_elec[ggid] += velecsum;
 325
 326         /* Increment number of inner iterations */
 327         inneriter                  += j_index_end - j_index_start;
 328
 329         /* Outer loop uses 31 flops */
 330     }
 331
 332     /* Increment number of outer iterations */
 333     outeriter        += nri;
 334
 335     /* Update outer/inner flops */
 336
 337     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_W3_VF,outeriter*31 + inneriter*96);
 338 }
 339 /*
 340  * Gromacs nonbonded kernel:   nb_kernel_ElecRFCut_VdwNone_GeomW3P1_F_c
 341  * Electrostatics interaction: ReactionField
 342  * VdW interaction:            None
 343  * Geometry:                   Water3-Particle
 344  * Calculate force/pot:        Force
 345  */
 346 void
 347 nb_kernel_ElecRFCut_VdwNone_GeomW3P1_F_c
 348                     (t_nblist                    * gmx_restrict       nlist,
 349                      rvec                        * gmx_restrict          xx,
 350                      rvec                        * gmx_restrict          ff,
 351                      t_forcerec                  * gmx_restrict          fr,
 352                      t_mdatoms                   * gmx_restrict     mdatoms,
 353                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
 354                      t_nrnb                      * gmx_restrict        nrnb)
 355 {
 356     int              i_shift_offset,i_coord_offset,j_coord_offset;
 357     int              j_index_start,j_index_end;
 358     int              nri,inr,ggid,iidx,jidx,jnr,outeriter,inneriter;
 359     real             shX,shY,shZ,tx,ty,tz,fscal,rcutoff,rcutoff2;
 360     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
 361     real             *shiftvec,*fshift,*x,*f;
 362     int              vdwioffset0;
 363     real             ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
 364     int              vdwioffset1;
 365     real             ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
 366     int              vdwioffset2;
 367     real             ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
 368     int              vdwjidx0;
 369     real             jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
 370     real             dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00,cexp1_00,cexp2_00;
 371     real             dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10,cexp1_10,cexp2_10;
 372     real             dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20,cexp1_20,cexp2_20;
 373     real             velec,felec,velecsum,facel,crf,krf,krf2;
 374     real             *charge;
 375
 376     x                = xx[0];
 377     f                = ff[0];
 378
 379     nri              = nlist->nri;
 380     iinr             = nlist->iinr;
 381     jindex           = nlist->jindex;
 382     jjnr             = nlist->jjnr;
 383     shiftidx         = nlist->shift;
 384     gid              = nlist->gid;
 385     shiftvec         = fr->shift_vec[0];
 386     fshift           = fr->fshift[0];
 387     facel            = fr->epsfac;
 388     charge           = mdatoms->chargeA;
 389     krf              = fr->ic->k_rf;
 390     krf2             = krf*2.0;
 391     crf              = fr->ic->c_rf;
 392
 393     /* Setup water-specific parameters */
 394     inr              = nlist->iinr[0];
 395     iq0              = facel*charge[inr+0];
 396     iq1              = facel*charge[inr+1];
 397     iq2              = facel*charge[inr+2];
 398
 399     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
 400     rcutoff          = fr->rcoulomb;
 401     rcutoff2         = rcutoff*rcutoff;
 402
 403     outeriter        = 0;
 404     inneriter        = 0;
 405
 406     /* Start outer loop over neighborlists */
 407     for(iidx=0; iidx<nri; iidx++)
 408     {
 409         /* Load shift vector for this list */
 410         i_shift_offset   = DIM*shiftidx[iidx];
 411         shX              = shiftvec[i_shift_offset+XX];
 412         shY              = shiftvec[i_shift_offset+YY];
 413         shZ              = shiftvec[i_shift_offset+ZZ];
 414
 415         /* Load limits for loop over neighbors */
 416         j_index_start    = jindex[iidx];
 417         j_index_end      = jindex[iidx+1];
 418
 419         /* Get outer coordinate index */
 420         inr              = iinr[iidx];
 421         i_coord_offset   = DIM*inr;
 422
 423         /* Load i particle coords and add shift vector */
 424         ix0              = shX + x[i_coord_offset+DIM*0+XX];
 425         iy0              = shY + x[i_coord_offset+DIM*0+YY];
 426         iz0              = shZ + x[i_coord_offset+DIM*0+ZZ];
 427         ix1              = shX + x[i_coord_offset+DIM*1+XX];
 428         iy1              = shY + x[i_coord_offset+DIM*1+YY];
 429         iz1              = shZ + x[i_coord_offset+DIM*1+ZZ];
 430         ix2              = shX + x[i_coord_offset+DIM*2+XX];
 431         iy2              = shY + x[i_coord_offset+DIM*2+YY];
 432         iz2              = shZ + x[i_coord_offset+DIM*2+ZZ];
 433
 434         fix0             = 0.0;
 435         fiy0             = 0.0;
 436         fiz0             = 0.0;
 437         fix1             = 0.0;
 438         fiy1             = 0.0;
 439         fiz1             = 0.0;
 440         fix2             = 0.0;
 441         fiy2             = 0.0;
 442         fiz2             = 0.0;
 443
 444         /* Start inner kernel loop */
 445         for(jidx=j_index_start; jidx<j_index_end; jidx++)
 446         {
 447             /* Get j neighbor index, and coordinate index */
 448             jnr              = jjnr[jidx];
 449             j_coord_offset   = DIM*jnr;
 450
 451             /* load j atom coordinates */
 452             jx0              = x[j_coord_offset+DIM*0+XX];
 453             jy0              = x[j_coord_offset+DIM*0+YY];
 454             jz0              = x[j_coord_offset+DIM*0+ZZ];
 455
 456             /* Calculate displacement vector */
 457             dx00             = ix0 - jx0;
 458             dy00             = iy0 - jy0;
 459             dz00             = iz0 - jz0;
 460             dx10             = ix1 - jx0;
 461             dy10             = iy1 - jy0;
 462             dz10             = iz1 - jz0;
 463             dx20             = ix2 - jx0;
 464             dy20             = iy2 - jy0;
 465             dz20             = iz2 - jz0;
 466
 467             /* Calculate squared distance and things based on it */
 468             rsq00            = dx00*dx00+dy00*dy00+dz00*dz00;
 469             rsq10            = dx10*dx10+dy10*dy10+dz10*dz10;
 470             rsq20            = dx20*dx20+dy20*dy20+dz20*dz20;
 471
 472             rinv00           = gmx_invsqrt(rsq00);
 473             rinv10           = gmx_invsqrt(rsq10);
 474             rinv20           = gmx_invsqrt(rsq20);
 475
 476             rinvsq00         = rinv00*rinv00;
 477             rinvsq10         = rinv10*rinv10;
 478             rinvsq20         = rinv20*rinv20;
 479
 480             /* Load parameters for j particles */
 481             jq0              = charge[jnr+0];
 482
 483             /**************************
 484              * CALCULATE INTERACTIONS *
 485              **************************/
 486
 487             if (rsq00<rcutoff2)
 488             {
 489
 490             qq00             = iq0*jq0;
 491
 492             /* REACTION-FIELD ELECTROSTATICS */
 493             felec            = qq00*(rinv00*rinvsq00-krf2);
 494
 495             fscal            = felec;
 496
 497             /* Calculate temporary vectorial force */
 498             tx               = fscal*dx00;
 499             ty               = fscal*dy00;
 500             tz               = fscal*dz00;
 501
 502             /* Update vectorial force */
 503             fix0            += tx;
 504             fiy0            += ty;
 505             fiz0            += tz;
 506             f[j_coord_offset+DIM*0+XX] -= tx;
 507             f[j_coord_offset+DIM*0+YY] -= ty;
 508             f[j_coord_offset+DIM*0+ZZ] -= tz;
 509
 510             }
 511
 512             /**************************
 513              * CALCULATE INTERACTIONS *
 514              **************************/
 515
 516             if (rsq10<rcutoff2)
 517             {
 518
 519             qq10             = iq1*jq0;
 520
 521             /* REACTION-FIELD ELECTROSTATICS */
 522             felec            = qq10*(rinv10*rinvsq10-krf2);
 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             qq20             = iq2*jq0;
 549
 550             /* REACTION-FIELD ELECTROSTATICS */
 551             felec            = qq20*(rinv20*rinvsq20-krf2);
 552
 553             fscal            = felec;
 554
 555             /* Calculate temporary vectorial force */
 556             tx               = fscal*dx20;
 557             ty               = fscal*dy20;
 558             tz               = fscal*dz20;
 559
 560             /* Update vectorial force */
 561             fix2            += tx;
 562             fiy2            += ty;
 563             fiz2            += tz;
 564             f[j_coord_offset+DIM*0+XX] -= tx;
 565             f[j_coord_offset+DIM*0+YY] -= ty;
 566             f[j_coord_offset+DIM*0+ZZ] -= tz;
 567
 568             }
 569
 570             /* Inner loop uses 81 flops */
 571         }
 572         /* End of innermost loop */
 573
 574         tx = ty = tz = 0;
 575         f[i_coord_offset+DIM*0+XX] += fix0;
 576         f[i_coord_offset+DIM*0+YY] += fiy0;
 577         f[i_coord_offset+DIM*0+ZZ] += fiz0;
 578         tx                         += fix0;
 579         ty                         += fiy0;
 580         tz                         += fiz0;
 581         f[i_coord_offset+DIM*1+XX] += fix1;
 582         f[i_coord_offset+DIM*1+YY] += fiy1;
 583         f[i_coord_offset+DIM*1+ZZ] += fiz1;
 584         tx                         += fix1;
 585         ty                         += fiy1;
 586         tz                         += fiz1;
 587         f[i_coord_offset+DIM*2+XX] += fix2;
 588         f[i_coord_offset+DIM*2+YY] += fiy2;
 589         f[i_coord_offset+DIM*2+ZZ] += fiz2;
 590         tx                         += fix2;
 591         ty                         += fiy2;
 592         tz                         += fiz2;
 593         fshift[i_shift_offset+XX]  += tx;
 594         fshift[i_shift_offset+YY]  += ty;
 595         fshift[i_shift_offset+ZZ]  += tz;
 596
 597         /* Increment number of inner iterations */
 598         inneriter                  += j_index_end - j_index_start;
 599
 600         /* Outer loop uses 30 flops */
 601     }
 602
 603     /* Increment number of outer iterations */
 604     outeriter        += nri;
 605
 606     /* Update outer/inner flops */
 607
 608     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_W3_F,outeriter*30 + inneriter*81);
 609 }