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