src/gmxlib/nonbonded/nb_kernel_c/nb_kernel_ElecRF_VdwBham_GeomW3W3_c.c

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