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