src/gmxlib/nonbonded/nb_kernel_power6/pwr6kernel334.F

   1 C
   2 C                This source code is part of
   3 C
   4 C                 G   R   O   M   A   C   S
   5 C
   6 C Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   7 C Copyright (c) 2001-2009, The GROMACS Development Team
   8 C
   9 C Gromacs is a library for molecular simulation and trajectory analysis,
  10 C written by Erik Lindahl, David van der Spoel, Berk Hess, and others - for
  11 C a full list of developers and information, check out http://www.gromacs.org
  12 C
  13 C This program is free software; you can redistribute it and/or modify it under
  14 C the terms of the GNU Lesser General Public License as published by the Free
  15 C Software Foundation; either version 2 of the License, or (at your option) any
  16 C later version.
  17 C As a special exception, you may use this file as part of a free software
  18 C library without restriction.  Specifically, if other files instantiate
  19 C templates or use macros or inline functions from this file, or you compile
  20 C this file and link it with other files to produce an executable, this
  21 C file does not by itself cause the resulting executable to be covered by
  22 C the GNU Lesser General Public License.
  23 C
  24 C In plain-speak: do not worry about classes/macros/templates either - only
  25 C changes to the library have to be LGPL, not an application linking with it.
  26 C
  27 C To help fund GROMACS development, we humbly ask that you cite
  28 C the papers people have written on it - you can find them on the website!
  29 C
  30
  31 #ifdef HAVE_CONFIG_H
  32 #  include<config.h>
  33 #endif
  34
  35 #ifdef GMX_DOUBLE
  36 #  define gmxreal real*8
  37 #else
  38 #  define gmxreal real*4
  39 #endif
  40
  41
  42
  43 C
  44 C Gromacs nonbonded kernel pwr6kernel334
  45 C Coulomb interaction:     Tabulated
  46 C VdW interaction:         Tabulated
  47 C water optimization:      pairs of TIP4P interactions
  48 C Calculate forces:        yes
  49 C
  50       subroutine pwr6kernel334(
  51      &                          nri,
  52      &                          iinr,
  53      &                          jindex,
  54      &                          jjnr,
  55      &                          shift,
  56      &                          shiftvec,
  57      &                          fshift,
  58      &                          gid,
  59      &                          pos,
  60      &                          faction,
  61      &                          charge,
  62      &                          facel,
  63      &                          krf,
  64      &                          crf,
  65      &                          Vc,
  66      &                          type,
  67      &                          ntype,
  68      &                          vdwparam,
  69      &                          Vvdw,
  70      &                          tabscale,
  71      &                          VFtab,
  72      &                          invsqrta,
  73      &                          dvda,
  74      &                          gbtabscale,
  75      &                          GBtab,
  76      &                          nthreads,
  77      &                          count,
  78      &                          mtx,
  79      &                          outeriter,
  80      &                          inneriter,
  81      &                          work)
  82       implicit      none
  83       integer*4     nri,iinr(*),jindex(*),jjnr(*),shift(*)
  84       gmxreal       shiftvec(*),fshift(*),pos(*),faction(*)
  85       integer*4     gid(*),type(*),ntype
  86       gmxreal       charge(*),facel,krf,crf,Vc(*),vdwparam(*)
  87       gmxreal       Vvdw(*),tabscale,VFtab(*)
  88       gmxreal       invsqrta(*),dvda(*),gbtabscale,GBtab(*)
  89       integer*4     nthreads,count,mtx,outeriter,inneriter
  90       gmxreal       work(*)
  91
  92       integer*4     n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid
  93       integer*4     nn0,nn1,nouter,ninner
  94       gmxreal       shX,shY,shZ
  95       gmxreal       fscal,tx,ty,tz
  96       gmxreal       qq,vcoul,vctot
  97       integer*4     tj
  98       gmxreal       Vvdw6,Vvdwtot
  99       gmxreal       Vvdw12
 100       gmxreal       r,rt,eps,eps2
 101       integer*4     n0,nnn
 102       gmxreal       Y,F,Geps,Heps2,Fp,VV
 103       gmxreal       FF
 104       gmxreal       fijC
 105       gmxreal       fijD,fijR
 106       gmxreal       ix1,iy1,iz1,fix1,fiy1,fiz1
 107       gmxreal       ix2,iy2,iz2,fix2,fiy2,fiz2
 108       gmxreal       ix3,iy3,iz3,fix3,fiy3,fiz3
 109       gmxreal       ix4,iy4,iz4,fix4,fiy4,fiz4
 110       gmxreal       jx1,jy1,jz1
 111       gmxreal       jx2,jy2,jz2,fjx2,fjy2,fjz2
 112       gmxreal       jx3,jy3,jz3,fjx3,fjy3,fjz3
 113       gmxreal       jx4,jy4,jz4,fjx4,fjy4,fjz4
 114       gmxreal       dx11,dy11,dz11,rsq11,rinv11
 115       gmxreal       dx22,dy22,dz22,rsq22,rinv22
 116       gmxreal       dx23,dy23,dz23,rsq23,rinv23
 117       gmxreal       dx24,dy24,dz24,rsq24,rinv24
 118       gmxreal       dx32,dy32,dz32,rsq32,rinv32
 119       gmxreal       dx33,dy33,dz33,rsq33,rinv33
 120       gmxreal       dx34,dy34,dz34,rsq34,rinv34
 121       gmxreal       dx42,dy42,dz42,rsq42,rinv42
 122       gmxreal       dx43,dy43,dz43,rsq43,rinv43
 123       gmxreal       dx44,dy44,dz44,rsq44,rinv44
 124       gmxreal       qH,qM,qqMM,qqMH,qqHH
 125       gmxreal       c6,c12
 126
 127
 128 C    Initialize water data
 129       ii               = iinr(1)+1
 130       qH               = charge(ii+1)
 131       qM               = charge(ii+3)
 132       qqMM             = facel*qM*qM
 133       qqMH             = facel*qM*qH
 134       qqHH             = facel*qH*qH
 135       tj               = 2*(ntype+1)*type(ii)+1
 136       c6               = vdwparam(tj)
 137       c12              = vdwparam(tj+1)
 138
 139
 140 C    Reset outer and inner iteration counters
 141       nouter           = 0
 142       ninner           = 0
 143
 144 C    Loop over thread workunits
 145    10 call pwr6kernelsync(mtx,count,nri,nthreads,nn0,nn1)
 146         if(nn1.gt.nri) nn1=nri
 147
 148 C      Start outer loop over neighborlists
 149
 150         do n=nn0+1,nn1
 151
 152 C        Load shift vector for this list
 153           is3              = 3*shift(n)+1
 154           shX              = shiftvec(is3)
 155           shY              = shiftvec(is3+1)
 156           shZ              = shiftvec(is3+2)
 157
 158 C        Load limits for loop over neighbors
 159           nj0              = jindex(n)+1
 160           nj1              = jindex(n+1)
 161
 162 C        Get outer coordinate index
 163           ii               = iinr(n)+1
 164           ii3              = 3*ii-2
 165
 166 C        Load i atom data, add shift vector
 167           ix1              = shX + pos(ii3+0)
 168           iy1              = shY + pos(ii3+1)
 169           iz1              = shZ + pos(ii3+2)
 170           ix2              = shX + pos(ii3+3)
 171           iy2              = shY + pos(ii3+4)
 172           iz2              = shZ + pos(ii3+5)
 173           ix3              = shX + pos(ii3+6)
 174           iy3              = shY + pos(ii3+7)
 175           iz3              = shZ + pos(ii3+8)
 176           ix4              = shX + pos(ii3+9)
 177           iy4              = shY + pos(ii3+10)
 178           iz4              = shZ + pos(ii3+11)
 179
 180 C        Zero the potential energy for this list
 181           vctot            = 0
 182           Vvdwtot          = 0
 183
 184 C        Clear i atom forces
 185           fix1             = 0
 186           fiy1             = 0
 187           fiz1             = 0
 188           fix2             = 0
 189           fiy2             = 0
 190           fiz2             = 0
 191           fix3             = 0
 192           fiy3             = 0
 193           fiz3             = 0
 194           fix4             = 0
 195           fiy4             = 0
 196           fiz4             = 0
 197
 198           do k=nj0,nj1
 199
 200 C          Get j neighbor index, and coordinate index
 201             jnr              = jjnr(k)+1
 202             j3               = 3*jnr-2
 203
 204 C          load j atom coordinates
 205             jx1              = pos(j3+0)
 206             jy1              = pos(j3+1)
 207             jz1              = pos(j3+2)
 208             jx2              = pos(j3+3)
 209             jy2              = pos(j3+4)
 210             jz2              = pos(j3+5)
 211             jx3              = pos(j3+6)
 212             jy3              = pos(j3+7)
 213             jz3              = pos(j3+8)
 214             jx4              = pos(j3+9)
 215             jy4              = pos(j3+10)
 216             jz4              = pos(j3+11)
 217
 218 C          Calculate distance
 219             dx11             = ix1 - jx1
 220             dy11             = iy1 - jy1
 221             dz11             = iz1 - jz1
 222             rsq11            = dx11*dx11+dy11*dy11+dz11*dz11
 223             dx22             = ix2 - jx2
 224             dy22             = iy2 - jy2
 225             dz22             = iz2 - jz2
 226             rsq22            = dx22*dx22+dy22*dy22+dz22*dz22
 227             dx23             = ix2 - jx3
 228             dy23             = iy2 - jy3
 229             dz23             = iz2 - jz3
 230             rsq23            = dx23*dx23+dy23*dy23+dz23*dz23
 231             dx24             = ix2 - jx4
 232             dy24             = iy2 - jy4
 233             dz24             = iz2 - jz4
 234             rsq24            = dx24*dx24+dy24*dy24+dz24*dz24
 235             dx32             = ix3 - jx2
 236             dy32             = iy3 - jy2
 237             dz32             = iz3 - jz2
 238             rsq32            = dx32*dx32+dy32*dy32+dz32*dz32
 239             dx33             = ix3 - jx3
 240             dy33             = iy3 - jy3
 241             dz33             = iz3 - jz3
 242             rsq33            = dx33*dx33+dy33*dy33+dz33*dz33
 243             dx34             = ix3 - jx4
 244             dy34             = iy3 - jy4
 245             dz34             = iz3 - jz4
 246             rsq34            = dx34*dx34+dy34*dy34+dz34*dz34
 247             dx42             = ix4 - jx2
 248             dy42             = iy4 - jy2
 249             dz42             = iz4 - jz2
 250             rsq42            = dx42*dx42+dy42*dy42+dz42*dz42
 251             dx43             = ix4 - jx3
 252             dy43             = iy4 - jy3
 253             dz43             = iz4 - jz3
 254             rsq43            = dx43*dx43+dy43*dy43+dz43*dz43
 255             dx44             = ix4 - jx4
 256             dy44             = iy4 - jy4
 257             dz44             = iz4 - jz4
 258             rsq44            = dx44*dx44+dy44*dy44+dz44*dz44
 259
 260 C          Calculate 1/r and 1/r2
 261
 262 C          PowerPC intrinsics 1/sqrt lookup table
 263 #ifndef GMX_BLUEGENE
 264             rinv11           = frsqrtes(rsq11)
 265 #else
 266             rinv11           = frsqrte(dble(rsq11))
 267 #endif
 268             rinv11           = (0.5*rinv11*(3.0-((rsq11*rinv11)
 269      &  *rinv11)))
 270 #ifdef GMX_DOUBLE
 271             rinv11           = (0.5*rinv11*(3.0-((rsq11*rinv11)
 272      &  *rinv11)))
 273 #endif
 274
 275 C          PowerPC intrinsics 1/sqrt lookup table
 276 #ifndef GMX_BLUEGENE
 277             rinv22           = frsqrtes(rsq22)
 278 #else
 279             rinv22           = frsqrte(dble(rsq22))
 280 #endif
 281             rinv22           = (0.5*rinv22*(3.0-((rsq22*rinv22)
 282      &  *rinv22)))
 283 #ifdef GMX_DOUBLE
 284             rinv22           = (0.5*rinv22*(3.0-((rsq22*rinv22)
 285      &  *rinv22)))
 286 #endif
 287
 288 C          PowerPC intrinsics 1/sqrt lookup table
 289 #ifndef GMX_BLUEGENE
 290             rinv23           = frsqrtes(rsq23)
 291 #else
 292             rinv23           = frsqrte(dble(rsq23))
 293 #endif
 294             rinv23           = (0.5*rinv23*(3.0-((rsq23*rinv23)
 295      &  *rinv23)))
 296 #ifdef GMX_DOUBLE
 297             rinv23           = (0.5*rinv23*(3.0-((rsq23*rinv23)
 298      &  *rinv23)))
 299 #endif
 300
 301 C          PowerPC intrinsics 1/sqrt lookup table
 302 #ifndef GMX_BLUEGENE
 303             rinv24           = frsqrtes(rsq24)
 304 #else
 305             rinv24           = frsqrte(dble(rsq24))
 306 #endif
 307             rinv24           = (0.5*rinv24*(3.0-((rsq24*rinv24)
 308      &  *rinv24)))
 309 #ifdef GMX_DOUBLE
 310             rinv24           = (0.5*rinv24*(3.0-((rsq24*rinv24)
 311      &  *rinv24)))
 312 #endif
 313
 314 C          PowerPC intrinsics 1/sqrt lookup table
 315 #ifndef GMX_BLUEGENE
 316             rinv32           = frsqrtes(rsq32)
 317 #else
 318             rinv32           = frsqrte(dble(rsq32))
 319 #endif
 320             rinv32           = (0.5*rinv32*(3.0-((rsq32*rinv32)
 321      &  *rinv32)))
 322 #ifdef GMX_DOUBLE
 323             rinv32           = (0.5*rinv32*(3.0-((rsq32*rinv32)
 324      &  *rinv32)))
 325 #endif
 326
 327 C          PowerPC intrinsics 1/sqrt lookup table
 328 #ifndef GMX_BLUEGENE
 329             rinv33           = frsqrtes(rsq33)
 330 #else
 331             rinv33           = frsqrte(dble(rsq33))
 332 #endif
 333             rinv33           = (0.5*rinv33*(3.0-((rsq33*rinv33)
 334      &  *rinv33)))
 335 #ifdef GMX_DOUBLE
 336             rinv33           = (0.5*rinv33*(3.0-((rsq33*rinv33)
 337      &  *rinv33)))
 338 #endif
 339
 340 C          PowerPC intrinsics 1/sqrt lookup table
 341 #ifndef GMX_BLUEGENE
 342             rinv34           = frsqrtes(rsq34)
 343 #else
 344             rinv34           = frsqrte(dble(rsq34))
 345 #endif
 346             rinv34           = (0.5*rinv34*(3.0-((rsq34*rinv34)
 347      &  *rinv34)))
 348 #ifdef GMX_DOUBLE
 349             rinv34           = (0.5*rinv34*(3.0-((rsq34*rinv34)
 350      &  *rinv34)))
 351 #endif
 352
 353 C          PowerPC intrinsics 1/sqrt lookup table
 354 #ifndef GMX_BLUEGENE
 355             rinv42           = frsqrtes(rsq42)
 356 #else
 357             rinv42           = frsqrte(dble(rsq42))
 358 #endif
 359             rinv42           = (0.5*rinv42*(3.0-((rsq42*rinv42)
 360      &  *rinv42)))
 361 #ifdef GMX_DOUBLE
 362             rinv42           = (0.5*rinv42*(3.0-((rsq42*rinv42)
 363      &  *rinv42)))
 364 #endif
 365
 366 C          PowerPC intrinsics 1/sqrt lookup table
 367 #ifndef GMX_BLUEGENE
 368             rinv43           = frsqrtes(rsq43)
 369 #else
 370             rinv43           = frsqrte(dble(rsq43))
 371 #endif
 372             rinv43           = (0.5*rinv43*(3.0-((rsq43*rinv43)
 373      &  *rinv43)))
 374 #ifdef GMX_DOUBLE
 375             rinv43           = (0.5*rinv43*(3.0-((rsq43*rinv43)
 376      &  *rinv43)))
 377 #endif
 378
 379 C          PowerPC intrinsics 1/sqrt lookup table
 380 #ifndef GMX_BLUEGENE
 381             rinv44           = frsqrtes(rsq44)
 382 #else
 383             rinv44           = frsqrte(dble(rsq44))
 384 #endif
 385             rinv44           = (0.5*rinv44*(3.0-((rsq44*rinv44)
 386      &  *rinv44)))
 387 #ifdef GMX_DOUBLE
 388             rinv44           = (0.5*rinv44*(3.0-((rsq44*rinv44)
 389      &  *rinv44)))
 390 #endif
 391
 392 C          Load parameters for j atom
 393
 394 C          Calculate table index
 395             r                = rsq11*rinv11
 396
 397 C          Calculate table index
 398             rt               = r*tabscale
 399             n0               = rt
 400             eps              = rt-n0
 401             eps2             = eps*eps
 402             nnn              = 12*n0+1
 403
 404 C          Tabulated VdW interaction - dispersion
 405             nnn              = nnn+4
 406             Y                = VFtab(nnn)
 407             F                = VFtab(nnn+1)
 408             Geps             = eps*VFtab(nnn+2)
 409             Heps2            = eps2*VFtab(nnn+3)
 410             Fp               = F+Geps+Heps2
 411             VV               = Y+eps*Fp
 412             FF               = Fp+Geps+2.0*Heps2
 413             Vvdw6            = c6*VV
 414             fijD             = c6*FF
 415
 416 C          Tabulated VdW interaction - repulsion
 417             nnn              = nnn+4
 418             Y                = VFtab(nnn)
 419             F                = VFtab(nnn+1)
 420             Geps             = eps*VFtab(nnn+2)
 421             Heps2            = eps2*VFtab(nnn+3)
 422             Fp               = F+Geps+Heps2
 423             VV               = Y+eps*Fp
 424             FF               = Fp+Geps+2.0*Heps2
 425             Vvdw12           = c12*VV
 426             fijR             = c12*FF
 427             Vvdwtot          = Vvdwtot+ Vvdw6 + Vvdw12
 428             fscal            = -((fijD+fijR)*tabscale)*rinv11
 429
 430 C          Calculate temporary vectorial force
 431             tx               = fscal*dx11
 432             ty               = fscal*dy11
 433             tz               = fscal*dz11
 434
 435 C          Increment i atom force
 436             fix1             = fix1 + tx
 437             fiy1             = fiy1 + ty
 438             fiz1             = fiz1 + tz
 439
 440 C          Decrement j atom force
 441             faction(j3+0)    = faction(j3+0) - tx
 442             faction(j3+1)    = faction(j3+1) - ty
 443             faction(j3+2)    = faction(j3+2) - tz
 444
 445 C          Load parameters for j atom
 446             qq               = qqHH
 447
 448 C          Calculate table index
 449             r                = rsq22*rinv22
 450
 451 C          Calculate table index
 452             rt               = r*tabscale
 453             n0               = rt
 454             eps              = rt-n0
 455             eps2             = eps*eps
 456             nnn              = 12*n0+1
 457
 458 C          Tabulated coulomb interaction
 459             Y                = VFtab(nnn)
 460             F                = VFtab(nnn+1)
 461             Geps             = eps*VFtab(nnn+2)
 462             Heps2            = eps2*VFtab(nnn+3)
 463             Fp               = F+Geps+Heps2
 464             VV               = Y+eps*Fp
 465             FF               = Fp+Geps+2.0*Heps2
 466             vcoul            = qq*VV
 467             fijC             = qq*FF
 468             vctot            = vctot + vcoul
 469             fscal            = -((fijC)*tabscale)*rinv22
 470
 471 C          Calculate temporary vectorial force
 472             tx               = fscal*dx22
 473             ty               = fscal*dy22
 474             tz               = fscal*dz22
 475
 476 C          Increment i atom force
 477             fix2             = fix2 + tx
 478             fiy2             = fiy2 + ty
 479             fiz2             = fiz2 + tz
 480
 481 C          Decrement j atom force
 482             fjx2             = faction(j3+3) - tx
 483             fjy2             = faction(j3+4) - ty
 484             fjz2             = faction(j3+5) - tz
 485
 486 C          Load parameters for j atom
 487             qq               = qqHH
 488
 489 C          Calculate table index
 490             r                = rsq23*rinv23
 491
 492 C          Calculate table index
 493             rt               = r*tabscale
 494             n0               = rt
 495             eps              = rt-n0
 496             eps2             = eps*eps
 497             nnn              = 12*n0+1
 498
 499 C          Tabulated coulomb interaction
 500             Y                = VFtab(nnn)
 501             F                = VFtab(nnn+1)
 502             Geps             = eps*VFtab(nnn+2)
 503             Heps2            = eps2*VFtab(nnn+3)
 504             Fp               = F+Geps+Heps2
 505             VV               = Y+eps*Fp
 506             FF               = Fp+Geps+2.0*Heps2
 507             vcoul            = qq*VV
 508             fijC             = qq*FF
 509             vctot            = vctot + vcoul
 510             fscal            = -((fijC)*tabscale)*rinv23
 511
 512 C          Calculate temporary vectorial force
 513             tx               = fscal*dx23
 514             ty               = fscal*dy23
 515             tz               = fscal*dz23
 516
 517 C          Increment i atom force
 518             fix2             = fix2 + tx
 519             fiy2             = fiy2 + ty
 520             fiz2             = fiz2 + tz
 521
 522 C          Decrement j atom force
 523             fjx3             = faction(j3+6) - tx
 524             fjy3             = faction(j3+7) - ty
 525             fjz3             = faction(j3+8) - tz
 526
 527 C          Load parameters for j atom
 528             qq               = qqMH
 529
 530 C          Calculate table index
 531             r                = rsq24*rinv24
 532
 533 C          Calculate table index
 534             rt               = r*tabscale
 535             n0               = rt
 536             eps              = rt-n0
 537             eps2             = eps*eps
 538             nnn              = 12*n0+1
 539
 540 C          Tabulated coulomb interaction
 541             Y                = VFtab(nnn)
 542             F                = VFtab(nnn+1)
 543             Geps             = eps*VFtab(nnn+2)
 544             Heps2            = eps2*VFtab(nnn+3)
 545             Fp               = F+Geps+Heps2
 546             VV               = Y+eps*Fp
 547             FF               = Fp+Geps+2.0*Heps2
 548             vcoul            = qq*VV
 549             fijC             = qq*FF
 550             vctot            = vctot + vcoul
 551             fscal            = -((fijC)*tabscale)*rinv24
 552
 553 C          Calculate temporary vectorial force
 554             tx               = fscal*dx24
 555             ty               = fscal*dy24
 556             tz               = fscal*dz24
 557
 558 C          Increment i atom force
 559             fix2             = fix2 + tx
 560             fiy2             = fiy2 + ty
 561             fiz2             = fiz2 + tz
 562
 563 C          Decrement j atom force
 564             fjx4             = faction(j3+9) - tx
 565             fjy4             = faction(j3+10) - ty
 566             fjz4             = faction(j3+11) - tz
 567
 568 C          Load parameters for j atom
 569             qq               = qqHH
 570
 571 C          Calculate table index
 572             r                = rsq32*rinv32
 573
 574 C          Calculate table index
 575             rt               = r*tabscale
 576             n0               = rt
 577             eps              = rt-n0
 578             eps2             = eps*eps
 579             nnn              = 12*n0+1
 580
 581 C          Tabulated coulomb interaction
 582             Y                = VFtab(nnn)
 583             F                = VFtab(nnn+1)
 584             Geps             = eps*VFtab(nnn+2)
 585             Heps2            = eps2*VFtab(nnn+3)
 586             Fp               = F+Geps+Heps2
 587             VV               = Y+eps*Fp
 588             FF               = Fp+Geps+2.0*Heps2
 589             vcoul            = qq*VV
 590             fijC             = qq*FF
 591             vctot            = vctot + vcoul
 592             fscal            = -((fijC)*tabscale)*rinv32
 593
 594 C          Calculate temporary vectorial force
 595             tx               = fscal*dx32
 596             ty               = fscal*dy32
 597             tz               = fscal*dz32
 598
 599 C          Increment i atom force
 600             fix3             = fix3 + tx
 601             fiy3             = fiy3 + ty
 602             fiz3             = fiz3 + tz
 603
 604 C          Decrement j atom force
 605             fjx2             = fjx2 - tx
 606             fjy2             = fjy2 - ty
 607             fjz2             = fjz2 - tz
 608
 609 C          Load parameters for j atom
 610             qq               = qqHH
 611
 612 C          Calculate table index
 613             r                = rsq33*rinv33
 614
 615 C          Calculate table index
 616             rt               = r*tabscale
 617             n0               = rt
 618             eps              = rt-n0
 619             eps2             = eps*eps
 620             nnn              = 12*n0+1
 621
 622 C          Tabulated coulomb interaction
 623             Y                = VFtab(nnn)
 624             F                = VFtab(nnn+1)
 625             Geps             = eps*VFtab(nnn+2)
 626             Heps2            = eps2*VFtab(nnn+3)
 627             Fp               = F+Geps+Heps2
 628             VV               = Y+eps*Fp
 629             FF               = Fp+Geps+2.0*Heps2
 630             vcoul            = qq*VV
 631             fijC             = qq*FF
 632             vctot            = vctot + vcoul
 633             fscal            = -((fijC)*tabscale)*rinv33
 634
 635 C          Calculate temporary vectorial force
 636             tx               = fscal*dx33
 637             ty               = fscal*dy33
 638             tz               = fscal*dz33
 639
 640 C          Increment i atom force
 641             fix3             = fix3 + tx
 642             fiy3             = fiy3 + ty
 643             fiz3             = fiz3 + tz
 644
 645 C          Decrement j atom force
 646             fjx3             = fjx3 - tx
 647             fjy3             = fjy3 - ty
 648             fjz3             = fjz3 - tz
 649
 650 C          Load parameters for j atom
 651             qq               = qqMH
 652
 653 C          Calculate table index
 654             r                = rsq34*rinv34
 655
 656 C          Calculate table index
 657             rt               = r*tabscale
 658             n0               = rt
 659             eps              = rt-n0
 660             eps2             = eps*eps
 661             nnn              = 12*n0+1
 662
 663 C          Tabulated coulomb interaction
 664             Y                = VFtab(nnn)
 665             F                = VFtab(nnn+1)
 666             Geps             = eps*VFtab(nnn+2)
 667             Heps2            = eps2*VFtab(nnn+3)
 668             Fp               = F+Geps+Heps2
 669             VV               = Y+eps*Fp
 670             FF               = Fp+Geps+2.0*Heps2
 671             vcoul            = qq*VV
 672             fijC             = qq*FF
 673             vctot            = vctot + vcoul
 674             fscal            = -((fijC)*tabscale)*rinv34
 675
 676 C          Calculate temporary vectorial force
 677             tx               = fscal*dx34
 678             ty               = fscal*dy34
 679             tz               = fscal*dz34
 680
 681 C          Increment i atom force
 682             fix3             = fix3 + tx
 683             fiy3             = fiy3 + ty
 684             fiz3             = fiz3 + tz
 685
 686 C          Decrement j atom force
 687             fjx4             = fjx4 - tx
 688             fjy4             = fjy4 - ty
 689             fjz4             = fjz4 - tz
 690
 691 C          Load parameters for j atom
 692             qq               = qqMH
 693
 694 C          Calculate table index
 695             r                = rsq42*rinv42
 696
 697 C          Calculate table index
 698             rt               = r*tabscale
 699             n0               = rt
 700             eps              = rt-n0
 701             eps2             = eps*eps
 702             nnn              = 12*n0+1
 703
 704 C          Tabulated coulomb interaction
 705             Y                = VFtab(nnn)
 706             F                = VFtab(nnn+1)
 707             Geps             = eps*VFtab(nnn+2)
 708             Heps2            = eps2*VFtab(nnn+3)
 709             Fp               = F+Geps+Heps2
 710             VV               = Y+eps*Fp
 711             FF               = Fp+Geps+2.0*Heps2
 712             vcoul            = qq*VV
 713             fijC             = qq*FF
 714             vctot            = vctot + vcoul
 715             fscal            = -((fijC)*tabscale)*rinv42
 716
 717 C          Calculate temporary vectorial force
 718             tx               = fscal*dx42
 719             ty               = fscal*dy42
 720             tz               = fscal*dz42
 721
 722 C          Increment i atom force
 723             fix4             = fix4 + tx
 724             fiy4             = fiy4 + ty
 725             fiz4             = fiz4 + tz
 726
 727 C          Decrement j atom force
 728             faction(j3+3)    = fjx2 - tx
 729             faction(j3+4)    = fjy2 - ty
 730             faction(j3+5)    = fjz2 - tz
 731
 732 C          Load parameters for j atom
 733             qq               = qqMH
 734
 735 C          Calculate table index
 736             r                = rsq43*rinv43
 737
 738 C          Calculate table index
 739             rt               = r*tabscale
 740             n0               = rt
 741             eps              = rt-n0
 742             eps2             = eps*eps
 743             nnn              = 12*n0+1
 744
 745 C          Tabulated coulomb interaction
 746             Y                = VFtab(nnn)
 747             F                = VFtab(nnn+1)
 748             Geps             = eps*VFtab(nnn+2)
 749             Heps2            = eps2*VFtab(nnn+3)
 750             Fp               = F+Geps+Heps2
 751             VV               = Y+eps*Fp
 752             FF               = Fp+Geps+2.0*Heps2
 753             vcoul            = qq*VV
 754             fijC             = qq*FF
 755             vctot            = vctot + vcoul
 756             fscal            = -((fijC)*tabscale)*rinv43
 757
 758 C          Calculate temporary vectorial force
 759             tx               = fscal*dx43
 760             ty               = fscal*dy43
 761             tz               = fscal*dz43
 762
 763 C          Increment i atom force
 764             fix4             = fix4 + tx
 765             fiy4             = fiy4 + ty
 766             fiz4             = fiz4 + tz
 767
 768 C          Decrement j atom force
 769             faction(j3+6)    = fjx3 - tx
 770             faction(j3+7)    = fjy3 - ty
 771             faction(j3+8)    = fjz3 - tz
 772
 773 C          Load parameters for j atom
 774             qq               = qqMM
 775
 776 C          Calculate table index
 777             r                = rsq44*rinv44
 778
 779 C          Calculate table index
 780             rt               = r*tabscale
 781             n0               = rt
 782             eps              = rt-n0
 783             eps2             = eps*eps
 784             nnn              = 12*n0+1
 785
 786 C          Tabulated coulomb interaction
 787             Y                = VFtab(nnn)
 788             F                = VFtab(nnn+1)
 789             Geps             = eps*VFtab(nnn+2)
 790             Heps2            = eps2*VFtab(nnn+3)
 791             Fp               = F+Geps+Heps2
 792             VV               = Y+eps*Fp
 793             FF               = Fp+Geps+2.0*Heps2
 794             vcoul            = qq*VV
 795             fijC             = qq*FF
 796             vctot            = vctot + vcoul
 797             fscal            = -((fijC)*tabscale)*rinv44
 798
 799 C          Calculate temporary vectorial force
 800             tx               = fscal*dx44
 801             ty               = fscal*dy44
 802             tz               = fscal*dz44
 803
 804 C          Increment i atom force
 805             fix4             = fix4 + tx
 806             fiy4             = fiy4 + ty
 807             fiz4             = fiz4 + tz
 808
 809 C          Decrement j atom force
 810             faction(j3+9)    = fjx4 - tx
 811             faction(j3+10)   = fjy4 - ty
 812             faction(j3+11)   = fjz4 - tz
 813
 814 C          Inner loop uses 433 flops/iteration
 815           end do
 816
 817
 818 C        Add i forces to mem and shifted force list
 819           faction(ii3+0)   = faction(ii3+0) + fix1
 820           faction(ii3+1)   = faction(ii3+1) + fiy1
 821           faction(ii3+2)   = faction(ii3+2) + fiz1
 822           faction(ii3+3)   = faction(ii3+3) + fix2
 823           faction(ii3+4)   = faction(ii3+4) + fiy2
 824           faction(ii3+5)   = faction(ii3+5) + fiz2
 825           faction(ii3+6)   = faction(ii3+6) + fix3
 826           faction(ii3+7)   = faction(ii3+7) + fiy3
 827           faction(ii3+8)   = faction(ii3+8) + fiz3
 828           faction(ii3+9)   = faction(ii3+9) + fix4
 829           faction(ii3+10)  = faction(ii3+10) + fiy4
 830           faction(ii3+11)  = faction(ii3+11) + fiz4
 831           fshift(is3)      = fshift(is3)+fix1+fix2+fix3+fix4
 832           fshift(is3+1)    = fshift(is3+1)+fiy1+fiy2+fiy3+fiy4
 833           fshift(is3+2)    = fshift(is3+2)+fiz1+fiz2+fiz3+fiz4
 834
 835 C        Add potential energies to the group for this list
 836           ggid             = gid(n)+1
 837           Vc(ggid)         = Vc(ggid) + vctot
 838           Vvdw(ggid)       = Vvdw(ggid) + Vvdwtot
 839
 840 C        Increment number of inner iterations
 841           ninner           = ninner + nj1 - nj0
 842
 843 C        Outer loop uses 38 flops/iteration
 844         end do
 845
 846
 847 C      Increment number of outer iterations
 848         nouter           = nouter + nn1 - nn0
 849       if(nn1.lt.nri) goto 10
 850
 851 C    Write outer/inner iteration count to pointers
 852       outeriter        = nouter
 853       inneriter        = ninner
 854       return
 855       end
 856
 857
 858
 859
 860
 861
 862 C
 863 C Gromacs nonbonded kernel pwr6kernel334nf
 864 C Coulomb interaction:     Tabulated
 865 C VdW interaction:         Tabulated
 866 C water optimization:      pairs of TIP4P interactions
 867 C Calculate forces:        no
 868 C
 869       subroutine pwr6kernel334nf(
 870      &                          nri,
 871      &                          iinr,
 872      &                          jindex,
 873      &                          jjnr,
 874      &                          shift,
 875      &                          shiftvec,
 876      &                          fshift,
 877      &                          gid,
 878      &                          pos,
 879      &                          faction,
 880      &                          charge,
 881      &                          facel,
 882      &                          krf,
 883      &                          crf,
 884      &                          Vc,
 885      &                          type,
 886      &                          ntype,
 887      &                          vdwparam,
 888      &                          Vvdw,
 889      &                          tabscale,
 890      &                          VFtab,
 891      &                          invsqrta,
 892      &                          dvda,
 893      &                          gbtabscale,
 894      &                          GBtab,
 895      &                          nthreads,
 896      &                          count,
 897      &                          mtx,
 898      &                          outeriter,
 899      &                          inneriter,
 900      &                          work)
 901       implicit      none
 902       integer*4     nri,iinr(*),jindex(*),jjnr(*),shift(*)
 903       gmxreal       shiftvec(*),fshift(*),pos(*),faction(*)
 904       integer*4     gid(*),type(*),ntype
 905       gmxreal       charge(*),facel,krf,crf,Vc(*),vdwparam(*)
 906       gmxreal       Vvdw(*),tabscale,VFtab(*)
 907       gmxreal       invsqrta(*),dvda(*),gbtabscale,GBtab(*)
 908       integer*4     nthreads,count,mtx,outeriter,inneriter
 909       gmxreal       work(*)
 910
 911       integer*4     n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid
 912       integer*4     nn0,nn1,nouter,ninner
 913       gmxreal       shX,shY,shZ
 914       gmxreal       qq,vcoul,vctot
 915       integer*4     tj
 916       gmxreal       Vvdw6,Vvdwtot
 917       gmxreal       Vvdw12
 918       gmxreal       r,rt,eps,eps2
 919       integer*4     n0,nnn
 920       gmxreal       Y,F,Geps,Heps2,Fp,VV
 921       gmxreal       ix1,iy1,iz1
 922       gmxreal       ix2,iy2,iz2
 923       gmxreal       ix3,iy3,iz3
 924       gmxreal       ix4,iy4,iz4
 925       gmxreal       jx1,jy1,jz1
 926       gmxreal       jx2,jy2,jz2
 927       gmxreal       jx3,jy3,jz3
 928       gmxreal       jx4,jy4,jz4
 929       gmxreal       dx11,dy11,dz11,rsq11,rinv11
 930       gmxreal       dx22,dy22,dz22,rsq22,rinv22
 931       gmxreal       dx23,dy23,dz23,rsq23,rinv23
 932       gmxreal       dx24,dy24,dz24,rsq24,rinv24
 933       gmxreal       dx32,dy32,dz32,rsq32,rinv32
 934       gmxreal       dx33,dy33,dz33,rsq33,rinv33
 935       gmxreal       dx34,dy34,dz34,rsq34,rinv34
 936       gmxreal       dx42,dy42,dz42,rsq42,rinv42
 937       gmxreal       dx43,dy43,dz43,rsq43,rinv43
 938       gmxreal       dx44,dy44,dz44,rsq44,rinv44
 939       gmxreal       qH,qM,qqMM,qqMH,qqHH
 940       gmxreal       c6,c12
 941
 942
 943 C    Initialize water data
 944       ii               = iinr(1)+1
 945       qH               = charge(ii+1)
 946       qM               = charge(ii+3)
 947       qqMM             = facel*qM*qM
 948       qqMH             = facel*qM*qH
 949       qqHH             = facel*qH*qH
 950       tj               = 2*(ntype+1)*type(ii)+1
 951       c6               = vdwparam(tj)
 952       c12              = vdwparam(tj+1)
 953
 954
 955 C    Reset outer and inner iteration counters
 956       nouter           = 0
 957       ninner           = 0
 958
 959 C    Loop over thread workunits
 960    10 call pwr6kernelsync(mtx,count,nri,nthreads,nn0,nn1)
 961         if(nn1.gt.nri) nn1=nri
 962
 963 C      Start outer loop over neighborlists
 964
 965         do n=nn0+1,nn1
 966
 967 C        Load shift vector for this list
 968           is3              = 3*shift(n)+1
 969           shX              = shiftvec(is3)
 970           shY              = shiftvec(is3+1)
 971           shZ              = shiftvec(is3+2)
 972
 973 C        Load limits for loop over neighbors
 974           nj0              = jindex(n)+1
 975           nj1              = jindex(n+1)
 976
 977 C        Get outer coordinate index
 978           ii               = iinr(n)+1
 979           ii3              = 3*ii-2
 980
 981 C        Load i atom data, add shift vector
 982           ix1              = shX + pos(ii3+0)
 983           iy1              = shY + pos(ii3+1)
 984           iz1              = shZ + pos(ii3+2)
 985           ix2              = shX + pos(ii3+3)
 986           iy2              = shY + pos(ii3+4)
 987           iz2              = shZ + pos(ii3+5)
 988           ix3              = shX + pos(ii3+6)
 989           iy3              = shY + pos(ii3+7)
 990           iz3              = shZ + pos(ii3+8)
 991           ix4              = shX + pos(ii3+9)
 992           iy4              = shY + pos(ii3+10)
 993           iz4              = shZ + pos(ii3+11)
 994
 995 C        Zero the potential energy for this list
 996           vctot            = 0
 997           Vvdwtot          = 0
 998
 999 C        Clear i atom forces
1000
1001           do k=nj0,nj1
1002
1003 C          Get j neighbor index, and coordinate index
1004             jnr              = jjnr(k)+1
1005             j3               = 3*jnr-2
1006
1007 C          load j atom coordinates
1008             jx1              = pos(j3+0)
1009             jy1              = pos(j3+1)
1010             jz1              = pos(j3+2)
1011             jx2              = pos(j3+3)
1012             jy2              = pos(j3+4)
1013             jz2              = pos(j3+5)
1014             jx3              = pos(j3+6)
1015             jy3              = pos(j3+7)
1016             jz3              = pos(j3+8)
1017             jx4              = pos(j3+9)
1018             jy4              = pos(j3+10)
1019             jz4              = pos(j3+11)
1020
1021 C          Calculate distance
1022             dx11             = ix1 - jx1
1023             dy11             = iy1 - jy1
1024             dz11             = iz1 - jz1
1025             rsq11            = dx11*dx11+dy11*dy11+dz11*dz11
1026             dx22             = ix2 - jx2
1027             dy22             = iy2 - jy2
1028             dz22             = iz2 - jz2
1029             rsq22            = dx22*dx22+dy22*dy22+dz22*dz22
1030             dx23             = ix2 - jx3
1031             dy23             = iy2 - jy3
1032             dz23             = iz2 - jz3
1033             rsq23            = dx23*dx23+dy23*dy23+dz23*dz23
1034             dx24             = ix2 - jx4
1035             dy24             = iy2 - jy4
1036             dz24             = iz2 - jz4
1037             rsq24            = dx24*dx24+dy24*dy24+dz24*dz24
1038             dx32             = ix3 - jx2
1039             dy32             = iy3 - jy2
1040             dz32             = iz3 - jz2
1041             rsq32            = dx32*dx32+dy32*dy32+dz32*dz32
1042             dx33             = ix3 - jx3
1043             dy33             = iy3 - jy3
1044             dz33             = iz3 - jz3
1045             rsq33            = dx33*dx33+dy33*dy33+dz33*dz33
1046             dx34             = ix3 - jx4
1047             dy34             = iy3 - jy4
1048             dz34             = iz3 - jz4
1049             rsq34            = dx34*dx34+dy34*dy34+dz34*dz34
1050             dx42             = ix4 - jx2
1051             dy42             = iy4 - jy2
1052             dz42             = iz4 - jz2
1053             rsq42            = dx42*dx42+dy42*dy42+dz42*dz42
1054             dx43             = ix4 - jx3
1055             dy43             = iy4 - jy3
1056             dz43             = iz4 - jz3
1057             rsq43            = dx43*dx43+dy43*dy43+dz43*dz43
1058             dx44             = ix4 - jx4
1059             dy44             = iy4 - jy4
1060             dz44             = iz4 - jz4
1061             rsq44            = dx44*dx44+dy44*dy44+dz44*dz44
1062
1063 C          Calculate 1/r and 1/r2
1064
1065 C          PowerPC intrinsics 1/sqrt lookup table
1066 #ifndef GMX_BLUEGENE
1067             rinv11           = frsqrtes(rsq11)
1068 #else
1069             rinv11           = frsqrte(dble(rsq11))
1070 #endif
1071             rinv11           = (0.5*rinv11*(3.0-((rsq11*rinv11)
1072      &  *rinv11)))
1073 #ifdef GMX_DOUBLE
1074             rinv11           = (0.5*rinv11*(3.0-((rsq11*rinv11)
1075      &  *rinv11)))
1076 #endif
1077
1078 C          PowerPC intrinsics 1/sqrt lookup table
1079 #ifndef GMX_BLUEGENE
1080             rinv22           = frsqrtes(rsq22)
1081 #else
1082             rinv22           = frsqrte(dble(rsq22))
1083 #endif
1084             rinv22           = (0.5*rinv22*(3.0-((rsq22*rinv22)
1085      &  *rinv22)))
1086 #ifdef GMX_DOUBLE
1087             rinv22           = (0.5*rinv22*(3.0-((rsq22*rinv22)
1088      &  *rinv22)))
1089 #endif
1090
1091 C          PowerPC intrinsics 1/sqrt lookup table
1092 #ifndef GMX_BLUEGENE
1093             rinv23           = frsqrtes(rsq23)
1094 #else
1095             rinv23           = frsqrte(dble(rsq23))
1096 #endif
1097             rinv23           = (0.5*rinv23*(3.0-((rsq23*rinv23)
1098      &  *rinv23)))
1099 #ifdef GMX_DOUBLE
1100             rinv23           = (0.5*rinv23*(3.0-((rsq23*rinv23)
1101      &  *rinv23)))
1102 #endif
1103
1104 C          PowerPC intrinsics 1/sqrt lookup table
1105 #ifndef GMX_BLUEGENE
1106             rinv24           = frsqrtes(rsq24)
1107 #else
1108             rinv24           = frsqrte(dble(rsq24))
1109 #endif
1110             rinv24           = (0.5*rinv24*(3.0-((rsq24*rinv24)
1111      &  *rinv24)))
1112 #ifdef GMX_DOUBLE
1113             rinv24           = (0.5*rinv24*(3.0-((rsq24*rinv24)
1114      &  *rinv24)))
1115 #endif
1116
1117 C          PowerPC intrinsics 1/sqrt lookup table
1118 #ifndef GMX_BLUEGENE
1119             rinv32           = frsqrtes(rsq32)
1120 #else
1121             rinv32           = frsqrte(dble(rsq32))
1122 #endif
1123             rinv32           = (0.5*rinv32*(3.0-((rsq32*rinv32)
1124      &  *rinv32)))
1125 #ifdef GMX_DOUBLE
1126             rinv32           = (0.5*rinv32*(3.0-((rsq32*rinv32)
1127      &  *rinv32)))
1128 #endif
1129
1130 C          PowerPC intrinsics 1/sqrt lookup table
1131 #ifndef GMX_BLUEGENE
1132             rinv33           = frsqrtes(rsq33)
1133 #else
1134             rinv33           = frsqrte(dble(rsq33))
1135 #endif
1136             rinv33           = (0.5*rinv33*(3.0-((rsq33*rinv33)
1137      &  *rinv33)))
1138 #ifdef GMX_DOUBLE
1139             rinv33           = (0.5*rinv33*(3.0-((rsq33*rinv33)
1140      &  *rinv33)))
1141 #endif
1142
1143 C          PowerPC intrinsics 1/sqrt lookup table
1144 #ifndef GMX_BLUEGENE
1145             rinv34           = frsqrtes(rsq34)
1146 #else
1147             rinv34           = frsqrte(dble(rsq34))
1148 #endif
1149             rinv34           = (0.5*rinv34*(3.0-((rsq34*rinv34)
1150      &  *rinv34)))
1151 #ifdef GMX_DOUBLE
1152             rinv34           = (0.5*rinv34*(3.0-((rsq34*rinv34)
1153      &  *rinv34)))
1154 #endif
1155
1156 C          PowerPC intrinsics 1/sqrt lookup table
1157 #ifndef GMX_BLUEGENE
1158             rinv42           = frsqrtes(rsq42)
1159 #else
1160             rinv42           = frsqrte(dble(rsq42))
1161 #endif
1162             rinv42           = (0.5*rinv42*(3.0-((rsq42*rinv42)
1163      &  *rinv42)))
1164 #ifdef GMX_DOUBLE
1165             rinv42           = (0.5*rinv42*(3.0-((rsq42*rinv42)
1166      &  *rinv42)))
1167 #endif
1168
1169 C          PowerPC intrinsics 1/sqrt lookup table
1170 #ifndef GMX_BLUEGENE
1171             rinv43           = frsqrtes(rsq43)
1172 #else
1173             rinv43           = frsqrte(dble(rsq43))
1174 #endif
1175             rinv43           = (0.5*rinv43*(3.0-((rsq43*rinv43)
1176      &  *rinv43)))
1177 #ifdef GMX_DOUBLE
1178             rinv43           = (0.5*rinv43*(3.0-((rsq43*rinv43)
1179      &  *rinv43)))
1180 #endif
1181
1182 C          PowerPC intrinsics 1/sqrt lookup table
1183 #ifndef GMX_BLUEGENE
1184             rinv44           = frsqrtes(rsq44)
1185 #else
1186             rinv44           = frsqrte(dble(rsq44))
1187 #endif
1188             rinv44           = (0.5*rinv44*(3.0-((rsq44*rinv44)
1189      &  *rinv44)))
1190 #ifdef GMX_DOUBLE
1191             rinv44           = (0.5*rinv44*(3.0-((rsq44*rinv44)
1192      &  *rinv44)))
1193 #endif
1194
1195 C          Load parameters for j atom
1196
1197 C          Calculate table index
1198             r                = rsq11*rinv11
1199
1200 C          Calculate table index
1201             rt               = r*tabscale
1202             n0               = rt
1203             eps              = rt-n0
1204             eps2             = eps*eps
1205             nnn              = 12*n0+1
1206
1207 C          Tabulated VdW interaction - dispersion
1208             nnn              = nnn+4
1209             Y                = VFtab(nnn)
1210             F                = VFtab(nnn+1)
1211             Geps             = eps*VFtab(nnn+2)
1212             Heps2            = eps2*VFtab(nnn+3)
1213             Fp               = F+Geps+Heps2
1214             VV               = Y+eps*Fp
1215             Vvdw6            = c6*VV
1216
1217 C          Tabulated VdW interaction - repulsion
1218             nnn              = nnn+4
1219             Y                = VFtab(nnn)
1220             F                = VFtab(nnn+1)
1221             Geps             = eps*VFtab(nnn+2)
1222             Heps2            = eps2*VFtab(nnn+3)
1223             Fp               = F+Geps+Heps2
1224             VV               = Y+eps*Fp
1225             Vvdw12           = c12*VV
1226             Vvdwtot          = Vvdwtot+ Vvdw6 + Vvdw12
1227
1228 C          Load parameters for j atom
1229             qq               = qqHH
1230
1231 C          Calculate table index
1232             r                = rsq22*rinv22
1233
1234 C          Calculate table index
1235             rt               = r*tabscale
1236             n0               = rt
1237             eps              = rt-n0
1238             eps2             = eps*eps
1239             nnn              = 12*n0+1
1240
1241 C          Tabulated coulomb interaction
1242             Y                = VFtab(nnn)
1243             F                = VFtab(nnn+1)
1244             Geps             = eps*VFtab(nnn+2)
1245             Heps2            = eps2*VFtab(nnn+3)
1246             Fp               = F+Geps+Heps2
1247             VV               = Y+eps*Fp
1248             vcoul            = qq*VV
1249             vctot            = vctot + vcoul
1250
1251 C          Load parameters for j atom
1252             qq               = qqHH
1253
1254 C          Calculate table index
1255             r                = rsq23*rinv23
1256
1257 C          Calculate table index
1258             rt               = r*tabscale
1259             n0               = rt
1260             eps              = rt-n0
1261             eps2             = eps*eps
1262             nnn              = 12*n0+1
1263
1264 C          Tabulated coulomb interaction
1265             Y                = VFtab(nnn)
1266             F                = VFtab(nnn+1)
1267             Geps             = eps*VFtab(nnn+2)
1268             Heps2            = eps2*VFtab(nnn+3)
1269             Fp               = F+Geps+Heps2
1270             VV               = Y+eps*Fp
1271             vcoul            = qq*VV
1272             vctot            = vctot + vcoul
1273
1274 C          Load parameters for j atom
1275             qq               = qqMH
1276
1277 C          Calculate table index
1278             r                = rsq24*rinv24
1279
1280 C          Calculate table index
1281             rt               = r*tabscale
1282             n0               = rt
1283             eps              = rt-n0
1284             eps2             = eps*eps
1285             nnn              = 12*n0+1
1286
1287 C          Tabulated coulomb interaction
1288             Y                = VFtab(nnn)
1289             F                = VFtab(nnn+1)
1290             Geps             = eps*VFtab(nnn+2)
1291             Heps2            = eps2*VFtab(nnn+3)
1292             Fp               = F+Geps+Heps2
1293             VV               = Y+eps*Fp
1294             vcoul            = qq*VV
1295             vctot            = vctot + vcoul
1296
1297 C          Load parameters for j atom
1298             qq               = qqHH
1299
1300 C          Calculate table index
1301             r                = rsq32*rinv32
1302
1303 C          Calculate table index
1304             rt               = r*tabscale
1305             n0               = rt
1306             eps              = rt-n0
1307             eps2             = eps*eps
1308             nnn              = 12*n0+1
1309
1310 C          Tabulated coulomb interaction
1311             Y                = VFtab(nnn)
1312             F                = VFtab(nnn+1)
1313             Geps             = eps*VFtab(nnn+2)
1314             Heps2            = eps2*VFtab(nnn+3)
1315             Fp               = F+Geps+Heps2
1316             VV               = Y+eps*Fp
1317             vcoul            = qq*VV
1318             vctot            = vctot + vcoul
1319
1320 C          Load parameters for j atom
1321             qq               = qqHH
1322
1323 C          Calculate table index
1324             r                = rsq33*rinv33
1325
1326 C          Calculate table index
1327             rt               = r*tabscale
1328             n0               = rt
1329             eps              = rt-n0
1330             eps2             = eps*eps
1331             nnn              = 12*n0+1
1332
1333 C          Tabulated coulomb interaction
1334             Y                = VFtab(nnn)
1335             F                = VFtab(nnn+1)
1336             Geps             = eps*VFtab(nnn+2)
1337             Heps2            = eps2*VFtab(nnn+3)
1338             Fp               = F+Geps+Heps2
1339             VV               = Y+eps*Fp
1340             vcoul            = qq*VV
1341             vctot            = vctot + vcoul
1342
1343 C          Load parameters for j atom
1344             qq               = qqMH
1345
1346 C          Calculate table index
1347             r                = rsq34*rinv34
1348
1349 C          Calculate table index
1350             rt               = r*tabscale
1351             n0               = rt
1352             eps              = rt-n0
1353             eps2             = eps*eps
1354             nnn              = 12*n0+1
1355
1356 C          Tabulated coulomb interaction
1357             Y                = VFtab(nnn)
1358             F                = VFtab(nnn+1)
1359             Geps             = eps*VFtab(nnn+2)
1360             Heps2            = eps2*VFtab(nnn+3)
1361             Fp               = F+Geps+Heps2
1362             VV               = Y+eps*Fp
1363             vcoul            = qq*VV
1364             vctot            = vctot + vcoul
1365
1366 C          Load parameters for j atom
1367             qq               = qqMH
1368
1369 C          Calculate table index
1370             r                = rsq42*rinv42
1371
1372 C          Calculate table index
1373             rt               = r*tabscale
1374             n0               = rt
1375             eps              = rt-n0
1376             eps2             = eps*eps
1377             nnn              = 12*n0+1
1378
1379 C          Tabulated coulomb interaction
1380             Y                = VFtab(nnn)
1381             F                = VFtab(nnn+1)
1382             Geps             = eps*VFtab(nnn+2)
1383             Heps2            = eps2*VFtab(nnn+3)
1384             Fp               = F+Geps+Heps2
1385             VV               = Y+eps*Fp
1386             vcoul            = qq*VV
1387             vctot            = vctot + vcoul
1388
1389 C          Load parameters for j atom
1390             qq               = qqMH
1391
1392 C          Calculate table index
1393             r                = rsq43*rinv43
1394
1395 C          Calculate table index
1396             rt               = r*tabscale
1397             n0               = rt
1398             eps              = rt-n0
1399             eps2             = eps*eps
1400             nnn              = 12*n0+1
1401
1402 C          Tabulated coulomb interaction
1403             Y                = VFtab(nnn)
1404             F                = VFtab(nnn+1)
1405             Geps             = eps*VFtab(nnn+2)
1406             Heps2            = eps2*VFtab(nnn+3)
1407             Fp               = F+Geps+Heps2
1408             VV               = Y+eps*Fp
1409             vcoul            = qq*VV
1410             vctot            = vctot + vcoul
1411
1412 C          Load parameters for j atom
1413             qq               = qqMM
1414
1415 C          Calculate table index
1416             r                = rsq44*rinv44
1417
1418 C          Calculate table index
1419             rt               = r*tabscale
1420             n0               = rt
1421             eps              = rt-n0
1422             eps2             = eps*eps
1423             nnn              = 12*n0+1
1424
1425 C          Tabulated coulomb interaction
1426             Y                = VFtab(nnn)
1427             F                = VFtab(nnn+1)
1428             Geps             = eps*VFtab(nnn+2)
1429             Heps2            = eps2*VFtab(nnn+3)
1430             Fp               = F+Geps+Heps2
1431             VV               = Y+eps*Fp
1432             vcoul            = qq*VV
1433             vctot            = vctot + vcoul
1434
1435 C          Inner loop uses 268 flops/iteration
1436           end do
1437
1438
1439 C        Add i forces to mem and shifted force list
1440
1441 C        Add potential energies to the group for this list
1442           ggid             = gid(n)+1
1443           Vc(ggid)         = Vc(ggid) + vctot
1444           Vvdw(ggid)       = Vvdw(ggid) + Vvdwtot
1445
1446 C        Increment number of inner iterations
1447           ninner           = ninner + nj1 - nj0
1448
1449 C        Outer loop uses 14 flops/iteration
1450         end do
1451
1452
1453 C      Increment number of outer iterations
1454         nouter           = nouter + nn1 - nn0
1455       if(nn1.lt.nri) goto 10
1456
1457 C    Write outer/inner iteration count to pointers
1458       outeriter        = nouter
1459       inneriter        = ninner
1460       return
1461       end
1462
1463
1464