src/gromacs/gmxlib/nonbonded/nb_kernel_c/nb_kernel_ElecRF_VdwNone_GeomW3W3_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 "config.h"
  39
  40 #include <math.h>
  41
  42 #include "../nb_kernel.h"
  43 #include "types/simple.h"
  44 #include "gromacs/math/vec.h"
  45 #include "nrnb.h"
  46
  47 /*
  48  * Gromacs nonbonded kernel:   nb_kernel_ElecRF_VdwNone_GeomW3W3_VF_c
  49  * Electrostatics interaction: ReactionField
  50  * VdW interaction:            None
  51  * Geometry:                   Water3-Water3
  52  * Calculate force/pot:        PotentialAndForce
  53  */
  54 void
  55 nb_kernel_ElecRF_VdwNone_GeomW3W3_VF_c
  56                     (t_nblist                    * gmx_restrict       nlist,
  57                      rvec                        * gmx_restrict          xx,
  58                      rvec                        * gmx_restrict          ff,
  59                      t_forcerec                  * gmx_restrict          fr,
  60                      t_mdatoms                   * gmx_restrict     mdatoms,
  61                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
  62                      t_nrnb                      * gmx_restrict        nrnb)
  63 {
  64     int              i_shift_offset,i_coord_offset,j_coord_offset;
  65     int              j_index_start,j_index_end;
  66     int              nri,inr,ggid,iidx,jidx,jnr,outeriter,inneriter;
  67     real             shX,shY,shZ,tx,ty,tz,fscal,rcutoff,rcutoff2;
  68     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
  69     real             *shiftvec,*fshift,*x,*f;
  70     int              vdwioffset0;
  71     real             ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
  72     int              vdwioffset1;
  73     real             ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
  74     int              vdwioffset2;
  75     real             ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
  76     int              vdwjidx0;
  77     real             jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
  78     int              vdwjidx1;
  79     real             jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
  80     int              vdwjidx2;
  81     real             jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
  82     real             dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00,cexp1_00,cexp2_00;
  83     real             dx01,dy01,dz01,rsq01,rinv01,rinvsq01,r01,qq01,c6_01,c12_01,cexp1_01,cexp2_01;
  84     real             dx02,dy02,dz02,rsq02,rinv02,rinvsq02,r02,qq02,c6_02,c12_02,cexp1_02,cexp2_02;
  85     real             dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10,cexp1_10,cexp2_10;
  86     real             dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11,cexp1_11,cexp2_11;
  87     real             dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12,cexp1_12,cexp2_12;
  88     real             dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20,cexp1_20,cexp2_20;
  89     real             dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21,cexp1_21,cexp2_21;
  90     real             dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22,cexp1_22,cexp2_22;
  91     real             velec,felec,velecsum,facel,crf,krf,krf2;
  92     real             *charge;
  93
  94     x                = xx[0];
  95     f                = ff[0];
  96
  97     nri              = nlist->nri;
  98     iinr             = nlist->iinr;
  99     jindex           = nlist->jindex;
 100     jjnr             = nlist->jjnr;
 101     shiftidx         = nlist->shift;
 102     gid              = nlist->gid;
 103     shiftvec         = fr->shift_vec[0];
 104     fshift           = fr->fshift[0];
 105     facel            = fr->epsfac;
 106     charge           = mdatoms->chargeA;
 107     krf              = fr->ic->k_rf;
 108     krf2             = krf*2.0;
 109     crf              = fr->ic->c_rf;
 110
 111     /* Setup water-specific parameters */
 112     inr              = nlist->iinr[0];
 113     iq0              = facel*charge[inr+0];
 114     iq1              = facel*charge[inr+1];
 115     iq2              = facel*charge[inr+2];
 116
 117     jq0              = charge[inr+0];
 118     jq1              = charge[inr+1];
 119     jq2              = charge[inr+2];
 120     qq00             = iq0*jq0;
 121     qq01             = iq0*jq1;
 122     qq02             = iq0*jq2;
 123     qq10             = iq1*jq0;
 124     qq11             = iq1*jq1;
 125     qq12             = iq1*jq2;
 126     qq20             = iq2*jq0;
 127     qq21             = iq2*jq1;
 128     qq22             = iq2*jq2;
 129
 130     outeriter        = 0;
 131     inneriter        = 0;
 132
 133     /* Start outer loop over neighborlists */
 134     for(iidx=0; iidx<nri; iidx++)
 135     {
 136         /* Load shift vector for this list */
 137         i_shift_offset   = DIM*shiftidx[iidx];
 138         shX              = shiftvec[i_shift_offset+XX];
 139         shY              = shiftvec[i_shift_offset+YY];
 140         shZ              = shiftvec[i_shift_offset+ZZ];
 141
 142         /* Load limits for loop over neighbors */
 143         j_index_start    = jindex[iidx];
 144         j_index_end      = jindex[iidx+1];
 145
 146         /* Get outer coordinate index */
 147         inr              = iinr[iidx];
 148         i_coord_offset   = DIM*inr;
 149
 150         /* Load i particle coords and add shift vector */
 151         ix0              = shX + x[i_coord_offset+DIM*0+XX];
 152         iy0              = shY + x[i_coord_offset+DIM*0+YY];
 153         iz0              = shZ + x[i_coord_offset+DIM*0+ZZ];
 154         ix1              = shX + x[i_coord_offset+DIM*1+XX];
 155         iy1              = shY + x[i_coord_offset+DIM*1+YY];
 156         iz1              = shZ + x[i_coord_offset+DIM*1+ZZ];
 157         ix2              = shX + x[i_coord_offset+DIM*2+XX];
 158         iy2              = shY + x[i_coord_offset+DIM*2+YY];
 159         iz2              = shZ + x[i_coord_offset+DIM*2+ZZ];
 160
 161         fix0             = 0.0;
 162         fiy0             = 0.0;
 163         fiz0             = 0.0;
 164         fix1             = 0.0;
 165         fiy1             = 0.0;
 166         fiz1             = 0.0;
 167         fix2             = 0.0;
 168         fiy2             = 0.0;
 169         fiz2             = 0.0;
 170
 171         /* Reset potential sums */
 172         velecsum         = 0.0;
 173
 174         /* Start inner kernel loop */
 175         for(jidx=j_index_start; jidx<j_index_end; jidx++)
 176         {
 177             /* Get j neighbor index, and coordinate index */
 178             jnr              = jjnr[jidx];
 179             j_coord_offset   = DIM*jnr;
 180
 181             /* load j atom coordinates */
 182             jx0              = x[j_coord_offset+DIM*0+XX];
 183             jy0              = x[j_coord_offset+DIM*0+YY];
 184             jz0              = x[j_coord_offset+DIM*0+ZZ];
 185             jx1              = x[j_coord_offset+DIM*1+XX];
 186             jy1              = x[j_coord_offset+DIM*1+YY];
 187             jz1              = x[j_coord_offset+DIM*1+ZZ];
 188             jx2              = x[j_coord_offset+DIM*2+XX];
 189             jy2              = x[j_coord_offset+DIM*2+YY];
 190             jz2              = x[j_coord_offset+DIM*2+ZZ];
 191
 192             /* Calculate displacement vector */
 193             dx00             = ix0 - jx0;
 194             dy00             = iy0 - jy0;
 195             dz00             = iz0 - jz0;
 196             dx01             = ix0 - jx1;
 197             dy01             = iy0 - jy1;
 198             dz01             = iz0 - jz1;
 199             dx02             = ix0 - jx2;
 200             dy02             = iy0 - jy2;
 201             dz02             = iz0 - jz2;
 202             dx10             = ix1 - jx0;
 203             dy10             = iy1 - jy0;
 204             dz10             = iz1 - jz0;
 205             dx11             = ix1 - jx1;
 206             dy11             = iy1 - jy1;
 207             dz11             = iz1 - jz1;
 208             dx12             = ix1 - jx2;
 209             dy12             = iy1 - jy2;
 210             dz12             = iz1 - jz2;
 211             dx20             = ix2 - jx0;
 212             dy20             = iy2 - jy0;
 213             dz20             = iz2 - jz0;
 214             dx21             = ix2 - jx1;
 215             dy21             = iy2 - jy1;
 216             dz21             = iz2 - jz1;
 217             dx22             = ix2 - jx2;
 218             dy22             = iy2 - jy2;
 219             dz22             = iz2 - jz2;
 220
 221             /* Calculate squared distance and things based on it */
 222             rsq00            = dx00*dx00+dy00*dy00+dz00*dz00;
 223             rsq01            = dx01*dx01+dy01*dy01+dz01*dz01;
 224             rsq02            = dx02*dx02+dy02*dy02+dz02*dz02;
 225             rsq10            = dx10*dx10+dy10*dy10+dz10*dz10;
 226             rsq11            = dx11*dx11+dy11*dy11+dz11*dz11;
 227             rsq12            = dx12*dx12+dy12*dy12+dz12*dz12;
 228             rsq20            = dx20*dx20+dy20*dy20+dz20*dz20;
 229             rsq21            = dx21*dx21+dy21*dy21+dz21*dz21;
 230             rsq22            = dx22*dx22+dy22*dy22+dz22*dz22;
 231
 232             rinv00           = gmx_invsqrt(rsq00);
 233             rinv01           = gmx_invsqrt(rsq01);
 234             rinv02           = gmx_invsqrt(rsq02);
 235             rinv10           = gmx_invsqrt(rsq10);
 236             rinv11           = gmx_invsqrt(rsq11);
 237             rinv12           = gmx_invsqrt(rsq12);
 238             rinv20           = gmx_invsqrt(rsq20);
 239             rinv21           = gmx_invsqrt(rsq21);
 240             rinv22           = gmx_invsqrt(rsq22);
 241
 242             rinvsq00         = rinv00*rinv00;
 243             rinvsq01         = rinv01*rinv01;
 244             rinvsq02         = rinv02*rinv02;
 245             rinvsq10         = rinv10*rinv10;
 246             rinvsq11         = rinv11*rinv11;
 247             rinvsq12         = rinv12*rinv12;
 248             rinvsq20         = rinv20*rinv20;
 249             rinvsq21         = rinv21*rinv21;
 250             rinvsq22         = rinv22*rinv22;
 251
 252             /**************************
 253              * CALCULATE INTERACTIONS *
 254              **************************/
 255
 256             /* REACTION-FIELD ELECTROSTATICS */
 257             velec            = qq00*(rinv00+krf*rsq00-crf);
 258             felec            = qq00*(rinv00*rinvsq00-krf2);
 259
 260             /* Update potential sums from outer loop */
 261             velecsum        += velec;
 262
 263             fscal            = felec;
 264
 265             /* Calculate temporary vectorial force */
 266             tx               = fscal*dx00;
 267             ty               = fscal*dy00;
 268             tz               = fscal*dz00;
 269
 270             /* Update vectorial force */
 271             fix0            += tx;
 272             fiy0            += ty;
 273             fiz0            += tz;
 274             f[j_coord_offset+DIM*0+XX] -= tx;
 275             f[j_coord_offset+DIM*0+YY] -= ty;
 276             f[j_coord_offset+DIM*0+ZZ] -= tz;
 277
 278             /**************************
 279              * CALCULATE INTERACTIONS *
 280              **************************/
 281
 282             /* REACTION-FIELD ELECTROSTATICS */
 283             velec            = qq01*(rinv01+krf*rsq01-crf);
 284             felec            = qq01*(rinv01*rinvsq01-krf2);
 285
 286             /* Update potential sums from outer loop */
 287             velecsum        += velec;
 288
 289             fscal            = felec;
 290
 291             /* Calculate temporary vectorial force */
 292             tx               = fscal*dx01;
 293             ty               = fscal*dy01;
 294             tz               = fscal*dz01;
 295
 296             /* Update vectorial force */
 297             fix0            += tx;
 298             fiy0            += ty;
 299             fiz0            += tz;
 300             f[j_coord_offset+DIM*1+XX] -= tx;
 301             f[j_coord_offset+DIM*1+YY] -= ty;
 302             f[j_coord_offset+DIM*1+ZZ] -= tz;
 303
 304             /**************************
 305              * CALCULATE INTERACTIONS *
 306              **************************/
 307
 308             /* REACTION-FIELD ELECTROSTATICS */
 309             velec            = qq02*(rinv02+krf*rsq02-crf);
 310             felec            = qq02*(rinv02*rinvsq02-krf2);
 311
 312             /* Update potential sums from outer loop */
 313             velecsum        += velec;
 314
 315             fscal            = felec;
 316
 317             /* Calculate temporary vectorial force */
 318             tx               = fscal*dx02;
 319             ty               = fscal*dy02;
 320             tz               = fscal*dz02;
 321
 322             /* Update vectorial force */
 323             fix0            += tx;
 324             fiy0            += ty;
 325             fiz0            += tz;
 326             f[j_coord_offset+DIM*2+XX] -= tx;
 327             f[j_coord_offset+DIM*2+YY] -= ty;
 328             f[j_coord_offset+DIM*2+ZZ] -= tz;
 329
 330             /**************************
 331              * CALCULATE INTERACTIONS *
 332              **************************/
 333
 334             /* REACTION-FIELD ELECTROSTATICS */
 335             velec            = qq10*(rinv10+krf*rsq10-crf);
 336             felec            = qq10*(rinv10*rinvsq10-krf2);
 337
 338             /* Update potential sums from outer loop */
 339             velecsum        += velec;
 340
 341             fscal            = felec;
 342
 343             /* Calculate temporary vectorial force */
 344             tx               = fscal*dx10;
 345             ty               = fscal*dy10;
 346             tz               = fscal*dz10;
 347
 348             /* Update vectorial force */
 349             fix1            += tx;
 350             fiy1            += ty;
 351             fiz1            += tz;
 352             f[j_coord_offset+DIM*0+XX] -= tx;
 353             f[j_coord_offset+DIM*0+YY] -= ty;
 354             f[j_coord_offset+DIM*0+ZZ] -= tz;
 355
 356             /**************************
 357              * CALCULATE INTERACTIONS *
 358              **************************/
 359
 360             /* REACTION-FIELD ELECTROSTATICS */
 361             velec            = qq11*(rinv11+krf*rsq11-crf);
 362             felec            = qq11*(rinv11*rinvsq11-krf2);
 363
 364             /* Update potential sums from outer loop */
 365             velecsum        += velec;
 366
 367             fscal            = felec;
 368
 369             /* Calculate temporary vectorial force */
 370             tx               = fscal*dx11;
 371             ty               = fscal*dy11;
 372             tz               = fscal*dz11;
 373
 374             /* Update vectorial force */
 375             fix1            += tx;
 376             fiy1            += ty;
 377             fiz1            += tz;
 378             f[j_coord_offset+DIM*1+XX] -= tx;
 379             f[j_coord_offset+DIM*1+YY] -= ty;
 380             f[j_coord_offset+DIM*1+ZZ] -= tz;
 381
 382             /**************************
 383              * CALCULATE INTERACTIONS *
 384              **************************/
 385
 386             /* REACTION-FIELD ELECTROSTATICS */
 387             velec            = qq12*(rinv12+krf*rsq12-crf);
 388             felec            = qq12*(rinv12*rinvsq12-krf2);
 389
 390             /* Update potential sums from outer loop */
 391             velecsum        += velec;
 392
 393             fscal            = felec;
 394
 395             /* Calculate temporary vectorial force */
 396             tx               = fscal*dx12;
 397             ty               = fscal*dy12;
 398             tz               = fscal*dz12;
 399
 400             /* Update vectorial force */
 401             fix1            += tx;
 402             fiy1            += ty;
 403             fiz1            += tz;
 404             f[j_coord_offset+DIM*2+XX] -= tx;
 405             f[j_coord_offset+DIM*2+YY] -= ty;
 406             f[j_coord_offset+DIM*2+ZZ] -= tz;
 407
 408             /**************************
 409              * CALCULATE INTERACTIONS *
 410              **************************/
 411
 412             /* REACTION-FIELD ELECTROSTATICS */
 413             velec            = qq20*(rinv20+krf*rsq20-crf);
 414             felec            = qq20*(rinv20*rinvsq20-krf2);
 415
 416             /* Update potential sums from outer loop */
 417             velecsum        += velec;
 418
 419             fscal            = felec;
 420
 421             /* Calculate temporary vectorial force */
 422             tx               = fscal*dx20;
 423             ty               = fscal*dy20;
 424             tz               = fscal*dz20;
 425
 426             /* Update vectorial force */
 427             fix2            += tx;
 428             fiy2            += ty;
 429             fiz2            += tz;
 430             f[j_coord_offset+DIM*0+XX] -= tx;
 431             f[j_coord_offset+DIM*0+YY] -= ty;
 432             f[j_coord_offset+DIM*0+ZZ] -= tz;
 433
 434             /**************************
 435              * CALCULATE INTERACTIONS *
 436              **************************/
 437
 438             /* REACTION-FIELD ELECTROSTATICS */
 439             velec            = qq21*(rinv21+krf*rsq21-crf);
 440             felec            = qq21*(rinv21*rinvsq21-krf2);
 441
 442             /* Update potential sums from outer loop */
 443             velecsum        += velec;
 444
 445             fscal            = felec;
 446
 447             /* Calculate temporary vectorial force */
 448             tx               = fscal*dx21;
 449             ty               = fscal*dy21;
 450             tz               = fscal*dz21;
 451
 452             /* Update vectorial force */
 453             fix2            += tx;
 454             fiy2            += ty;
 455             fiz2            += tz;
 456             f[j_coord_offset+DIM*1+XX] -= tx;
 457             f[j_coord_offset+DIM*1+YY] -= ty;
 458             f[j_coord_offset+DIM*1+ZZ] -= tz;
 459
 460             /**************************
 461              * CALCULATE INTERACTIONS *
 462              **************************/
 463
 464             /* REACTION-FIELD ELECTROSTATICS */
 465             velec            = qq22*(rinv22+krf*rsq22-crf);
 466             felec            = qq22*(rinv22*rinvsq22-krf2);
 467
 468             /* Update potential sums from outer loop */
 469             velecsum        += velec;
 470
 471             fscal            = felec;
 472
 473             /* Calculate temporary vectorial force */
 474             tx               = fscal*dx22;
 475             ty               = fscal*dy22;
 476             tz               = fscal*dz22;
 477
 478             /* Update vectorial force */
 479             fix2            += tx;
 480             fiy2            += ty;
 481             fiz2            += tz;
 482             f[j_coord_offset+DIM*2+XX] -= tx;
 483             f[j_coord_offset+DIM*2+YY] -= ty;
 484             f[j_coord_offset+DIM*2+ZZ] -= tz;
 485
 486             /* Inner loop uses 279 flops */
 487         }
 488         /* End of innermost loop */
 489
 490         tx = ty = tz = 0;
 491         f[i_coord_offset+DIM*0+XX] += fix0;
 492         f[i_coord_offset+DIM*0+YY] += fiy0;
 493         f[i_coord_offset+DIM*0+ZZ] += fiz0;
 494         tx                         += fix0;
 495         ty                         += fiy0;
 496         tz                         += fiz0;
 497         f[i_coord_offset+DIM*1+XX] += fix1;
 498         f[i_coord_offset+DIM*1+YY] += fiy1;
 499         f[i_coord_offset+DIM*1+ZZ] += fiz1;
 500         tx                         += fix1;
 501         ty                         += fiy1;
 502         tz                         += fiz1;
 503         f[i_coord_offset+DIM*2+XX] += fix2;
 504         f[i_coord_offset+DIM*2+YY] += fiy2;
 505         f[i_coord_offset+DIM*2+ZZ] += fiz2;
 506         tx                         += fix2;
 507         ty                         += fiy2;
 508         tz                         += fiz2;
 509         fshift[i_shift_offset+XX]  += tx;
 510         fshift[i_shift_offset+YY]  += ty;
 511         fshift[i_shift_offset+ZZ]  += tz;
 512
 513         ggid                        = gid[iidx];
 514         /* Update potential energies */
 515         kernel_data->energygrp_elec[ggid] += velecsum;
 516
 517         /* Increment number of inner iterations */
 518         inneriter                  += j_index_end - j_index_start;
 519
 520         /* Outer loop uses 31 flops */
 521     }
 522
 523     /* Increment number of outer iterations */
 524     outeriter        += nri;
 525
 526     /* Update outer/inner flops */
 527
 528     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_W3W3_VF,outeriter*31 + inneriter*279);
 529 }
 530 /*
 531  * Gromacs nonbonded kernel:   nb_kernel_ElecRF_VdwNone_GeomW3W3_F_c
 532  * Electrostatics interaction: ReactionField
 533  * VdW interaction:            None
 534  * Geometry:                   Water3-Water3
 535  * Calculate force/pot:        Force
 536  */
 537 void
 538 nb_kernel_ElecRF_VdwNone_GeomW3W3_F_c
 539                     (t_nblist                    * gmx_restrict       nlist,
 540                      rvec                        * gmx_restrict          xx,
 541                      rvec                        * gmx_restrict          ff,
 542                      t_forcerec                  * gmx_restrict          fr,
 543                      t_mdatoms                   * gmx_restrict     mdatoms,
 544                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
 545                      t_nrnb                      * gmx_restrict        nrnb)
 546 {
 547     int              i_shift_offset,i_coord_offset,j_coord_offset;
 548     int              j_index_start,j_index_end;
 549     int              nri,inr,ggid,iidx,jidx,jnr,outeriter,inneriter;
 550     real             shX,shY,shZ,tx,ty,tz,fscal,rcutoff,rcutoff2;
 551     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
 552     real             *shiftvec,*fshift,*x,*f;
 553     int              vdwioffset0;
 554     real             ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
 555     int              vdwioffset1;
 556     real             ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
 557     int              vdwioffset2;
 558     real             ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
 559     int              vdwjidx0;
 560     real             jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
 561     int              vdwjidx1;
 562     real             jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
 563     int              vdwjidx2;
 564     real             jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
 565     real             dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00,cexp1_00,cexp2_00;
 566     real             dx01,dy01,dz01,rsq01,rinv01,rinvsq01,r01,qq01,c6_01,c12_01,cexp1_01,cexp2_01;
 567     real             dx02,dy02,dz02,rsq02,rinv02,rinvsq02,r02,qq02,c6_02,c12_02,cexp1_02,cexp2_02;
 568     real             dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10,cexp1_10,cexp2_10;
 569     real             dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11,cexp1_11,cexp2_11;
 570     real             dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12,cexp1_12,cexp2_12;
 571     real             dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20,cexp1_20,cexp2_20;
 572     real             dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21,cexp1_21,cexp2_21;
 573     real             dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22,cexp1_22,cexp2_22;
 574     real             velec,felec,velecsum,facel,crf,krf,krf2;
 575     real             *charge;
 576
 577     x                = xx[0];
 578     f                = ff[0];
 579
 580     nri              = nlist->nri;
 581     iinr             = nlist->iinr;
 582     jindex           = nlist->jindex;
 583     jjnr             = nlist->jjnr;
 584     shiftidx         = nlist->shift;
 585     gid              = nlist->gid;
 586     shiftvec         = fr->shift_vec[0];
 587     fshift           = fr->fshift[0];
 588     facel            = fr->epsfac;
 589     charge           = mdatoms->chargeA;
 590     krf              = fr->ic->k_rf;
 591     krf2             = krf*2.0;
 592     crf              = fr->ic->c_rf;
 593
 594     /* Setup water-specific parameters */
 595     inr              = nlist->iinr[0];
 596     iq0              = facel*charge[inr+0];
 597     iq1              = facel*charge[inr+1];
 598     iq2              = facel*charge[inr+2];
 599
 600     jq0              = charge[inr+0];
 601     jq1              = charge[inr+1];
 602     jq2              = charge[inr+2];
 603     qq00             = iq0*jq0;
 604     qq01             = iq0*jq1;
 605     qq02             = iq0*jq2;
 606     qq10             = iq1*jq0;
 607     qq11             = iq1*jq1;
 608     qq12             = iq1*jq2;
 609     qq20             = iq2*jq0;
 610     qq21             = iq2*jq1;
 611     qq22             = iq2*jq2;
 612
 613     outeriter        = 0;
 614     inneriter        = 0;
 615
 616     /* Start outer loop over neighborlists */
 617     for(iidx=0; iidx<nri; iidx++)
 618     {
 619         /* Load shift vector for this list */
 620         i_shift_offset   = DIM*shiftidx[iidx];
 621         shX              = shiftvec[i_shift_offset+XX];
 622         shY              = shiftvec[i_shift_offset+YY];
 623         shZ              = shiftvec[i_shift_offset+ZZ];
 624
 625         /* Load limits for loop over neighbors */
 626         j_index_start    = jindex[iidx];
 627         j_index_end      = jindex[iidx+1];
 628
 629         /* Get outer coordinate index */
 630         inr              = iinr[iidx];
 631         i_coord_offset   = DIM*inr;
 632
 633         /* Load i particle coords and add shift vector */
 634         ix0              = shX + x[i_coord_offset+DIM*0+XX];
 635         iy0              = shY + x[i_coord_offset+DIM*0+YY];
 636         iz0              = shZ + x[i_coord_offset+DIM*0+ZZ];
 637         ix1              = shX + x[i_coord_offset+DIM*1+XX];
 638         iy1              = shY + x[i_coord_offset+DIM*1+YY];
 639         iz1              = shZ + x[i_coord_offset+DIM*1+ZZ];
 640         ix2              = shX + x[i_coord_offset+DIM*2+XX];
 641         iy2              = shY + x[i_coord_offset+DIM*2+YY];
 642         iz2              = shZ + x[i_coord_offset+DIM*2+ZZ];
 643
 644         fix0             = 0.0;
 645         fiy0             = 0.0;
 646         fiz0             = 0.0;
 647         fix1             = 0.0;
 648         fiy1             = 0.0;
 649         fiz1             = 0.0;
 650         fix2             = 0.0;
 651         fiy2             = 0.0;
 652         fiz2             = 0.0;
 653
 654         /* Start inner kernel loop */
 655         for(jidx=j_index_start; jidx<j_index_end; jidx++)
 656         {
 657             /* Get j neighbor index, and coordinate index */
 658             jnr              = jjnr[jidx];
 659             j_coord_offset   = DIM*jnr;
 660
 661             /* load j atom coordinates */
 662             jx0              = x[j_coord_offset+DIM*0+XX];
 663             jy0              = x[j_coord_offset+DIM*0+YY];
 664             jz0              = x[j_coord_offset+DIM*0+ZZ];
 665             jx1              = x[j_coord_offset+DIM*1+XX];
 666             jy1              = x[j_coord_offset+DIM*1+YY];
 667             jz1              = x[j_coord_offset+DIM*1+ZZ];
 668             jx2              = x[j_coord_offset+DIM*2+XX];
 669             jy2              = x[j_coord_offset+DIM*2+YY];
 670             jz2              = x[j_coord_offset+DIM*2+ZZ];
 671
 672             /* Calculate displacement vector */
 673             dx00             = ix0 - jx0;
 674             dy00             = iy0 - jy0;
 675             dz00             = iz0 - jz0;
 676             dx01             = ix0 - jx1;
 677             dy01             = iy0 - jy1;
 678             dz01             = iz0 - jz1;
 679             dx02             = ix0 - jx2;
 680             dy02             = iy0 - jy2;
 681             dz02             = iz0 - jz2;
 682             dx10             = ix1 - jx0;
 683             dy10             = iy1 - jy0;
 684             dz10             = iz1 - jz0;
 685             dx11             = ix1 - jx1;
 686             dy11             = iy1 - jy1;
 687             dz11             = iz1 - jz1;
 688             dx12             = ix1 - jx2;
 689             dy12             = iy1 - jy2;
 690             dz12             = iz1 - jz2;
 691             dx20             = ix2 - jx0;
 692             dy20             = iy2 - jy0;
 693             dz20             = iz2 - jz0;
 694             dx21             = ix2 - jx1;
 695             dy21             = iy2 - jy1;
 696             dz21             = iz2 - jz1;
 697             dx22             = ix2 - jx2;
 698             dy22             = iy2 - jy2;
 699             dz22             = iz2 - jz2;
 700
 701             /* Calculate squared distance and things based on it */
 702             rsq00            = dx00*dx00+dy00*dy00+dz00*dz00;
 703             rsq01            = dx01*dx01+dy01*dy01+dz01*dz01;
 704             rsq02            = dx02*dx02+dy02*dy02+dz02*dz02;
 705             rsq10            = dx10*dx10+dy10*dy10+dz10*dz10;
 706             rsq11            = dx11*dx11+dy11*dy11+dz11*dz11;
 707             rsq12            = dx12*dx12+dy12*dy12+dz12*dz12;
 708             rsq20            = dx20*dx20+dy20*dy20+dz20*dz20;
 709             rsq21            = dx21*dx21+dy21*dy21+dz21*dz21;
 710             rsq22            = dx22*dx22+dy22*dy22+dz22*dz22;
 711
 712             rinv00           = gmx_invsqrt(rsq00);
 713             rinv01           = gmx_invsqrt(rsq01);
 714             rinv02           = gmx_invsqrt(rsq02);
 715             rinv10           = gmx_invsqrt(rsq10);
 716             rinv11           = gmx_invsqrt(rsq11);
 717             rinv12           = gmx_invsqrt(rsq12);
 718             rinv20           = gmx_invsqrt(rsq20);
 719             rinv21           = gmx_invsqrt(rsq21);
 720             rinv22           = gmx_invsqrt(rsq22);
 721
 722             rinvsq00         = rinv00*rinv00;
 723             rinvsq01         = rinv01*rinv01;
 724             rinvsq02         = rinv02*rinv02;
 725             rinvsq10         = rinv10*rinv10;
 726             rinvsq11         = rinv11*rinv11;
 727             rinvsq12         = rinv12*rinv12;
 728             rinvsq20         = rinv20*rinv20;
 729             rinvsq21         = rinv21*rinv21;
 730             rinvsq22         = rinv22*rinv22;
 731
 732             /**************************
 733              * CALCULATE INTERACTIONS *
 734              **************************/
 735
 736             /* REACTION-FIELD ELECTROSTATICS */
 737             felec            = qq00*(rinv00*rinvsq00-krf2);
 738
 739             fscal            = felec;
 740
 741             /* Calculate temporary vectorial force */
 742             tx               = fscal*dx00;
 743             ty               = fscal*dy00;
 744             tz               = fscal*dz00;
 745
 746             /* Update vectorial force */
 747             fix0            += tx;
 748             fiy0            += ty;
 749             fiz0            += tz;
 750             f[j_coord_offset+DIM*0+XX] -= tx;
 751             f[j_coord_offset+DIM*0+YY] -= ty;
 752             f[j_coord_offset+DIM*0+ZZ] -= tz;
 753
 754             /**************************
 755              * CALCULATE INTERACTIONS *
 756              **************************/
 757
 758             /* REACTION-FIELD ELECTROSTATICS */
 759             felec            = qq01*(rinv01*rinvsq01-krf2);
 760
 761             fscal            = felec;
 762
 763             /* Calculate temporary vectorial force */
 764             tx               = fscal*dx01;
 765             ty               = fscal*dy01;
 766             tz               = fscal*dz01;
 767
 768             /* Update vectorial force */
 769             fix0            += tx;
 770             fiy0            += ty;
 771             fiz0            += tz;
 772             f[j_coord_offset+DIM*1+XX] -= tx;
 773             f[j_coord_offset+DIM*1+YY] -= ty;
 774             f[j_coord_offset+DIM*1+ZZ] -= tz;
 775
 776             /**************************
 777              * CALCULATE INTERACTIONS *
 778              **************************/
 779
 780             /* REACTION-FIELD ELECTROSTATICS */
 781             felec            = qq02*(rinv02*rinvsq02-krf2);
 782
 783             fscal            = felec;
 784
 785             /* Calculate temporary vectorial force */
 786             tx               = fscal*dx02;
 787             ty               = fscal*dy02;
 788             tz               = fscal*dz02;
 789
 790             /* Update vectorial force */
 791             fix0            += tx;
 792             fiy0            += ty;
 793             fiz0            += tz;
 794             f[j_coord_offset+DIM*2+XX] -= tx;
 795             f[j_coord_offset+DIM*2+YY] -= ty;
 796             f[j_coord_offset+DIM*2+ZZ] -= tz;
 797
 798             /**************************
 799              * CALCULATE INTERACTIONS *
 800              **************************/
 801
 802             /* REACTION-FIELD ELECTROSTATICS */
 803             felec            = qq10*(rinv10*rinvsq10-krf2);
 804
 805             fscal            = felec;
 806
 807             /* Calculate temporary vectorial force */
 808             tx               = fscal*dx10;
 809             ty               = fscal*dy10;
 810             tz               = fscal*dz10;
 811
 812             /* Update vectorial force */
 813             fix1            += tx;
 814             fiy1            += ty;
 815             fiz1            += tz;
 816             f[j_coord_offset+DIM*0+XX] -= tx;
 817             f[j_coord_offset+DIM*0+YY] -= ty;
 818             f[j_coord_offset+DIM*0+ZZ] -= tz;
 819
 820             /**************************
 821              * CALCULATE INTERACTIONS *
 822              **************************/
 823
 824             /* REACTION-FIELD ELECTROSTATICS */
 825             felec            = qq11*(rinv11*rinvsq11-krf2);
 826
 827             fscal            = felec;
 828
 829             /* Calculate temporary vectorial force */
 830             tx               = fscal*dx11;
 831             ty               = fscal*dy11;
 832             tz               = fscal*dz11;
 833
 834             /* Update vectorial force */
 835             fix1            += tx;
 836             fiy1            += ty;
 837             fiz1            += tz;
 838             f[j_coord_offset+DIM*1+XX] -= tx;
 839             f[j_coord_offset+DIM*1+YY] -= ty;
 840             f[j_coord_offset+DIM*1+ZZ] -= tz;
 841
 842             /**************************
 843              * CALCULATE INTERACTIONS *
 844              **************************/
 845
 846             /* REACTION-FIELD ELECTROSTATICS */
 847             felec            = qq12*(rinv12*rinvsq12-krf2);
 848
 849             fscal            = felec;
 850
 851             /* Calculate temporary vectorial force */
 852             tx               = fscal*dx12;
 853             ty               = fscal*dy12;
 854             tz               = fscal*dz12;
 855
 856             /* Update vectorial force */
 857             fix1            += tx;
 858             fiy1            += ty;
 859             fiz1            += tz;
 860             f[j_coord_offset+DIM*2+XX] -= tx;
 861             f[j_coord_offset+DIM*2+YY] -= ty;
 862             f[j_coord_offset+DIM*2+ZZ] -= tz;
 863
 864             /**************************
 865              * CALCULATE INTERACTIONS *
 866              **************************/
 867
 868             /* REACTION-FIELD ELECTROSTATICS */
 869             felec            = qq20*(rinv20*rinvsq20-krf2);
 870
 871             fscal            = felec;
 872
 873             /* Calculate temporary vectorial force */
 874             tx               = fscal*dx20;
 875             ty               = fscal*dy20;
 876             tz               = fscal*dz20;
 877
 878             /* Update vectorial force */
 879             fix2            += tx;
 880             fiy2            += ty;
 881             fiz2            += tz;
 882             f[j_coord_offset+DIM*0+XX] -= tx;
 883             f[j_coord_offset+DIM*0+YY] -= ty;
 884             f[j_coord_offset+DIM*0+ZZ] -= tz;
 885
 886             /**************************
 887              * CALCULATE INTERACTIONS *
 888              **************************/
 889
 890             /* REACTION-FIELD ELECTROSTATICS */
 891             felec            = qq21*(rinv21*rinvsq21-krf2);
 892
 893             fscal            = felec;
 894
 895             /* Calculate temporary vectorial force */
 896             tx               = fscal*dx21;
 897             ty               = fscal*dy21;
 898             tz               = fscal*dz21;
 899
 900             /* Update vectorial force */
 901             fix2            += tx;
 902             fiy2            += ty;
 903             fiz2            += tz;
 904             f[j_coord_offset+DIM*1+XX] -= tx;
 905             f[j_coord_offset+DIM*1+YY] -= ty;
 906             f[j_coord_offset+DIM*1+ZZ] -= tz;
 907
 908             /**************************
 909              * CALCULATE INTERACTIONS *
 910              **************************/
 911
 912             /* REACTION-FIELD ELECTROSTATICS */
 913             felec            = qq22*(rinv22*rinvsq22-krf2);
 914
 915             fscal            = felec;
 916
 917             /* Calculate temporary vectorial force */
 918             tx               = fscal*dx22;
 919             ty               = fscal*dy22;
 920             tz               = fscal*dz22;
 921
 922             /* Update vectorial force */
 923             fix2            += tx;
 924             fiy2            += ty;
 925             fiz2            += tz;
 926             f[j_coord_offset+DIM*2+XX] -= tx;
 927             f[j_coord_offset+DIM*2+YY] -= ty;
 928             f[j_coord_offset+DIM*2+ZZ] -= tz;
 929
 930             /* Inner loop uses 234 flops */
 931         }
 932         /* End of innermost loop */
 933
 934         tx = ty = tz = 0;
 935         f[i_coord_offset+DIM*0+XX] += fix0;
 936         f[i_coord_offset+DIM*0+YY] += fiy0;
 937         f[i_coord_offset+DIM*0+ZZ] += fiz0;
 938         tx                         += fix0;
 939         ty                         += fiy0;
 940         tz                         += fiz0;
 941         f[i_coord_offset+DIM*1+XX] += fix1;
 942         f[i_coord_offset+DIM*1+YY] += fiy1;
 943         f[i_coord_offset+DIM*1+ZZ] += fiz1;
 944         tx                         += fix1;
 945         ty                         += fiy1;
 946         tz                         += fiz1;
 947         f[i_coord_offset+DIM*2+XX] += fix2;
 948         f[i_coord_offset+DIM*2+YY] += fiy2;
 949         f[i_coord_offset+DIM*2+ZZ] += fiz2;
 950         tx                         += fix2;
 951         ty                         += fiy2;
 952         tz                         += fiz2;
 953         fshift[i_shift_offset+XX]  += tx;
 954         fshift[i_shift_offset+YY]  += ty;
 955         fshift[i_shift_offset+ZZ]  += tz;
 956
 957         /* Increment number of inner iterations */
 958         inneriter                  += j_index_end - j_index_start;
 959
 960         /* Outer loop uses 30 flops */
 961     }
 962
 963     /* Increment number of outer iterations */
 964     outeriter        += nri;
 965
 966     /* Update outer/inner flops */
 967
 968     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_W3W3_F,outeriter*30 + inneriter*234);
 969 }