src/gromacs/gmxlib/nonbonded/nb_kernel_c/nb_kernel_ElecRF_VdwLJ_GeomW4W4_c.c

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