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